problem with trsm and dtrsm

[Don_Ritchie]

I'm using Visual Fortran Composer XE for Windows on a machinerunning Windows 7 64-bit with i7-860 quad core and 8GB RAM.

In running some timings on the blas routines trsm and/or dtrsm, I've encountered a problem in getting correct results for moderate sized matrices. Everything looks good for matrices smaller than ca. 100x100, but when I try, for example, a 500x500 matrix, the 'maximum error' is ~1E+13.

Here's the code: <>

[bash]program quick_checks
use blas95

    REAL*8, ALLOCATABLE :: A(:,:),B(:,:),C(:,:),U(:,:)
    real*8 temp,t1,temp2
    character ans
    
    ans='y'
    do while((ans.eq.'y').or.(ans.eq.'Y'))
        print *,"Matrix size?"
        read *,N
    
        ALLOCATE(A(N,N),B(N,N),C(N,N),U(N,N))
        do i=1,N
            U(i,i)=1.0
            A(i,i)=1.0
            do j=i+1,N
               call random_number(A(i,j))
               call random_number(temp) !make A unit upper triangular matrix
               if (temp.gt.0.5) then
                   A(i,j)=-A(i,j)
               end if
               A(j,i)=0.0
               U(i,j)=0.0 
               U(j,i)=0.0
            end do
        end do
        temp=1.0
        t1=second()
        call trsm(A,U) !on output, U should be inverse of A
!        call dtrsm('L','U','N','U',N,N,temp,A,N,U,N)
        t1=second()-t1
    
        print *,"trsm execution time=",t1
        print *

        C=matmul(U,A)   !should be identity

        temp=0.0
        do i=1,N
            do j=1,N
                temp2=C(i,j)
                if (i.eq.j) then
                    temp2=temp2-1.0
                end if
                if (abs(temp2).gt.temp) then
                    temp=abs(temp2)
                    ibig=i
                    jbig=j
                end if
            end do
        end do

        print *,"maximum error=",temp," at i=",ibig," j=",jbig

        deallocate(B,C,A,U)
    
        print *,"Continue?"
        read *,ans
    end do
    
end program[/bash]

Either I'm doing something wrong, or there's a bug in these routines. I'd appreciate any help.

TIA

[Don_Ritchie]

I'd very much appreciate knowing if anyone has checked this.

[mecej4 · Accepted solution]

What is wrong is not what you are doing or the routines, but your expectations of such a calculation. You are obtaining the inverse of a large random matrix. Unless you can claim that the matrix has a small condition number, you should expect large errors when you compare A.inv(A) with I.

If all you are interested in is the timing of the MKL routines, your code has served its purpose. Do not try to extract more information from the exercise. Note, however, that pivoting is affected by matrix contents, so that the speed may vary slightly with the numbers in the matrix.

If you really want to know why you got large errors, use the MKL routines and evaluate the condition numbers of your random matrices and examine the errors vis-a-vis the condition numbers.

[Don_Ritchie]

Thank you. Lesson learned!! Now, I have tolearn how to evaluate "condition numbers".