Landed up in more disbelieving stuff. The way i was testing that subroutine (A) was printing the matrix(argument) that was passed on from B and then used that in the test program by reading. Now, it turns out that in the same program if i do this, that is, just print the matrix to a file and then immediately read, it works fine. something hard to believe but happened. Could it be that while writing and reading process (used only standard format i.e *) the digits to the order of 10^-5 or less could have changed but due to the sensitivity of the matrix caused the issue? any similar experience? pasting the subroutine below (i dont feel now as if there is a problem in it though yesterday said so), which basically uses a lapack subroutine. also below that the sample matrix tested which works fine is given for a rough idea of the magnitude of the terms. a portion of B which now works(?) is also pasted. Thanks for the quick replies.
any thoughts or suggestions?
is there a better way to share the code than pasting here. new to this forum.
--------------------
subroutine dginv(A,M,N,minmn,LWORK,Ainv)
IMPLICIT NONE
integer M,N,minmn,LWORK
real *8 eps
parameter(eps=1.0d-12)
real *8 A(M,N)
integer LDA,INFO,I,j,k,LDU
real *8 L(M,N),U(M,M),S(minmn),VT(N,N),V(N,N),WORK(LWORK)
integer LDVT
real *8 VSin(N,M),Ainv(N,M),tol,Sinv(minmn),Smax,UT(M,M)
character *1 JOBU,JOBVT
integer debug
! U(LDU,M),LDVT >= N,VT(LDVT,N)
! S(min(M,N)) work(LWORK)
! LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
JOBU='A'
JOBVT='A'
LDVT=N
LDA=MAX(M,1)
LDU=M
debug=1
write(*,*) '--------dginv Routine starts here-------------'
write(*,*) 'Inputs to the routine: '
write(*,*) M,N,minmn,LWORK
write(*,*) 'Matrix received is:'
do i=1,M
write(*,'(1X,100G14.5)') (A(i,j),j=1,N)
enddo
call DGESVD( JOBU, JOBVT, M, N, A, LDA, S,
; U, LDU, VT, LDVT,WORK, LWORK, INFO )
write(*,*) ' INFO : ',INFO
Smax=0.0d0
do i=1,minmn
if(abs(S(i))>Smax) Smax=S(i)
enddo
tol = eps*max(M,N)*Smax
do i=1,minmn
if(dabs(S(i))>tol) then
Sinv(i)=1./S(i)
else
Sinv(i)=0.0d0
endif
enddo
if(debug.eq.1) then
write(*,*) "Sinv:"
do i=1,minmn
write(*,*) Sinv(i)
enddo
endif
if(debug.eq.1) write(*,*) "V:"
do i=1,N
do j=1,N
V(i,j)=VT(j,i)
enddo
if(debug.eq.1) write(*,'(1X,100G14.5)') (V(i,j),j=1,N)
enddo
if(debug.eq.1) write(*,*) "UT:"
do i=1,M
do j=1,M
UT(i,j)=U(j,i)
enddo
if(debug.eq.1) write(*,'(1X,100G14.5)') (UT(i,j),j=1,M)
enddo
if(debug.eq.1) write(*,*) "Vsin:"
do i=1,N
do j=1,M
if(j<=minmn) then
Vsin(i,j)=V(i,j)*Sinv(j)
else
Vsin(i,j)=0.0d0
endif
enddo
if(debug.eq.1) write(*,'(1X,100G14.5)') (Vsin(i,j),j=1,M)
enddo
do i=1,N
do j=1,M
Ainv(i,j)=0.0d0
do k=1,M
Ainv(i,j)=Ainv(i,j)+Vsin(i,k)*UT(k,j)
enddo
enddo
enddo
write(*,*) 'Matrix leaving is:'
do i=1,N
write(*,'(1X,100G14.5)') (Ainv(i,j),j=1,M)
enddo
write(*,*) '--------dginv Routine ends here-------------'
! stop
end
--------------------
0.24835E-01 0.22708E-01 0.23572E-01 0.24288E-01 0.23835E-01 0.21405E-01 0.23744E-01 0.22450E-01
0.23406E-01 0.21617E-01 0.22436E-01 0.23198E-01 0.23406E-01 0.21274E-01 0.23596E-01 0.22342E-01
0.22977E-01 0.21527E-01 0.22300E-01 0.23107E-01 0.22977E-01 0.21144E-01 0.23448E-01 0.22235E-01
0.22548E-01 0.21436E-01 0.22164E-01 0.23015E-01 0.22548E-01 0.21014E-01 0.23298E-01 0.22127E-01
-0.33707E-02 -0.26382E-02 -0.39155E-02 -0.26141E-02 -0.23707E-02 -0.27780E-02 -0.32584E-02 -0.20972E-02
-0.34510E-02 -0.25935E-02 -0.38830E-02 -0.25915E-02 -0.34510E-02 -0.37291E-02 -0.42424E-02 -0.30630E-02
-0.35313E-02 -0.25904E-02 -0.38915E-02 -0.26107E-02 -0.35313E-02 -0.37212E-02 -0.42672E-02 -0.30702E-02
-0.36117E-02 -0.26317E-02 -0.39435E-02 -0.26743E-02 -0.36117E-02 -0.37569E-02 -0.43354E-02 -0.31215E-02
------------------------
code snippet for B calling A
if(debug==1) then
open(88,file="inputMat.txt")
write(*,*) '-----------------------Pmat----------------------'
do i=1,rows
write(*,'(1X,100G14.5)') (Pmat(i,j),j=1,cols)
write(88,'(1X,100G14.5)') (Pmat(i,j),j=1,cols)
enddo
write(*,*) '-----------------------Pmat----------------------'
write(*,*) 'Pmat -pinv descr: ',rows,cols,minmn,LWORK
close(88)
endif
! call dginv(Pmat,rows,cols,minmn,LWORK,Pinv)
if(debug==1) then
write(*,*) 'trial calculation'
open(99,file="inputMat.txt",iostat=IO)
do i=1,rows
read(99,*) (Pmat(i,j),j=1,cols)
enddo
write(*,*) 'PMAT read:'
do i=1,rows
write(*,'(1X,100G14.5)') (Pmat(i,j),j=1,cols)
enddo
close(99)
endif
call dginv(Pmat,rows,cols,minmn,LWORK,Pinv)