fortran - LAPACK inversion routine strangely mixes up all variables -



fortran - LAPACK inversion routine strangely mixes up all variables -

i'm programming in fortran , trying utilize dgetri matrix inverter lapack package:

http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html

but strangely seems messing variables. in simple example, matrix initialised @ origin of programme changes dgetri applied, though dgetri doesn't involve a…

can tell me going on? thanks!

program solvelinear real(8), dimension(2,2) :: a,ainv real(8),allocatable :: work(:) integer :: info,lwork,i integer,dimension(2) :: ipiv info=0 lwork=10000 allocate(work(lwork)) work=0 ipiv=0 = reshape((/1,-1,3,3/),(/2,2/)) ainv = reshape((/1,-1,3,3/),(/2,2/)) phone call dgetri(2,ainv,2,ipiv,work,lwork,info) print*,"a = " i=1,2 print*,a(i,:) end end programme solve linear

this output:

= 1.0000000000000000 0.0000000000000000 -1.0000000000000000 0.33333333333333331

you have compute lu decomposition before calling dgetri. using double version of routines lu has computed dgetrf (zgetrf complex version).

the next code compiles , returns right values

program solvelinear implicit none real(8), dimension(3,3) :: a,ainv,m,lu real(8),allocatable :: work(:) integer :: info,lwork integer,dimension(3) :: ipiv integer :: info=0 lwork=100 allocate(work(lwork)) work=0 ipiv=0 = reshape((/1,-1,3,3,1,1,1,1,1,1/),(/3,3/)) !-- lu factorisation lu = phone call dgetrf(3,3,lu,3,ipiv,info) !-- inversion of matrix using lu ainv=lu phone call dgetri(3,ainv,3,ipiv,work,lwork,info) !-- computation of a^-1 * check inverse m = matmul(ainv,a) print*,"m = " i=1,3 print*,m(i,:) end end programme solvelinear

output

m = 1.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000 -5.5511151231257827e-017 -4.4408920985006262e-016 -2.2204460492503131e-016 1.0000000000000000

by way, original code returns expected unchanged value of when compiled gfortran 4.8.2

fortran lapack linear-equation

Comments

Post a Comment

Popular posts from this blog

Delphi change the assembly code of a running process -

Delphi change the assembly code of a running process - i'm trying alter address 00741fa5 has push test.009e721c . alter push test.009e71c8 . procedure callback; asm force $9e71c8 end; procedure tform2.btn1click(sender: tobject); var ppid :dword; pprocess : cardinal; begin getwindowthreadprocessid(findwindow(nil,pchar('test')), @ppid); pprocess := openprocess(process_all_access, false, ppid); injectasm(pprocess, $741fa5, 10, integer(@callback), 0); end; here injectasm function (with error checking) function injectasm(process: longword; injectaddress, injectsize, codeaddress, codesize : integer): pointer; var csize, rsize : integer; replaced : array of byte; jmp : array [0..4] of byte; jmpaddress : integer; nopv : byte; : integer; nbr: ulong_ptr; begin //injectsize must equal or greater 5, because need space our //(far) jmp wwxxyyzz instruction //if there's no space, inject in place of few instructions if injectsize ...
Read more

json - Hibernate and Jackson (java.lang.IllegalStateException: Cannot call sendError() after the response has been committed) -

json - Hibernate and Jackson (java.lang.IllegalStateException: Cannot call sendError() after the response has been committed) - hello have problem hibernate , jackson. have 2 pojos @entity @table(name = "users") public class user implements serializable { static final long serialversionuid = 4235637833262722256l; @id @generatedvalue private int id; @column(name = "creation_date") private date creationdate; @column(name = "first_name") private string firstname; @column(name = "last_name") private string lastname; @version @column(name = "version") private int version; @column(name = "e_mail") private string email; @column(name = "enabled") private boolean enabled; @column(name = "login") private string login; @column(name = "password") private string password; @jsonmanagedreference("role") @onetomany(mappedby = "user", cascade = cascadetype.all,fetch = fetchtyp...
Read more

C++ 11 "class" keyword -

C++ 11 "class" keyword - i've started using c++ 11 little bit more , have few questions on special uses of class keyword. know it's used declare class, there's 2 instances have seen don't understand: method<class t>(); and class class_name *var; why precede typename keyword class in first example, , keyword pointer in sec example? that known elaborated type specifier , necessary when class name "shadowed" or "hidden" , need explicit. class t { }; // love of god don't t t; t t2; if compiler smart, give these warnings: main.cpp:15:5: error: must utilize 'class' tag refer type 't' in scope t t2; ^ class main.cpp:14:7: note: class 't' hidden non-type declaration of 't' here t t; ^ if class not defined, deed forwards declaration. for example: template <typename t> void method() { using type = typename t::value_type; } method...
Read more