Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations Mike Lewis on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

subroutine giving different results

Status
Not open for further replies.

ranjith1009

Programmer
Aug 13, 2009
4
US
I have a subroutine let us call it A(x,y,z).
Case 1: A is called from another subroutine say B (B is just one of the subroutines the main function uses-so a lot more .f files are compiled together along with it). The results i get are incorrect.
Case 2: But when i call the same subroutine A(x,y,z) with same arguments from a test program i get correct results.

There are no "common" variables accessed by A. The compiler used is f95. have tried using option -g but with same result.
The simplest check possible tried is: printing all the inputs and outputs of A(x,y,z) inside A. all inputs match in both cases (1 and 2) but output does not.
Will appreciate any help in solving this.
 
Do you have something like
Code:
program main
   integer x
   call fred(1)
   x = 1
   print *, x
end program

subroutine fred (xx)
   integer xx
   xx = 22
   return
end subroutine
 
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)
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top