topic for symmetric linear systems in Intel® Fortran Compiler

Cannot invert a matrix with IMSL routine LINDS

Pratap__Sangeeta — Mon, 29 Feb 2016 22:38:15 GMT

I am trying to use a standard IMSL routine LINDS to invert a positive definite matrix,.

I installed the IMSL library based on these instructions.

https://software.intel.com/en-us/articles/installing-and-using-the-imsl-libraries

Then I included the lines

INCLUDE 'link_fnl_shared.h'

use linds_int

in my program, before setting up the matrix to be inverted

The program compiles without error, but when I debug it it, it exits without inverting the matrix. It does not trigger a break point, but it just stops.

Any help will be greatly appreciated. Thank you very much

This would likely happen if

TimP — Mon, 29 Feb 2016 23:05:21 GMT

This would likely happen if you didn't set PATH to cover all dll you need at run time. If your IMSL is the one bought as an option with ifort, that should happen if you open the ifort cmd window. dependencywalker may be informative.

No, I am using Intel Visual

Pratap__Sangeeta — Tue, 01 Mar 2016 00:38:18 GMT

No, I am using Intel Visual Fortran for Windows. Not sure where to look..

Thanks!

When you run a Fortran

Steven_L_Intel1 — Tue, 01 Mar 2016 01:13:42 GMT

When you run a Fortran program in the Visual Studio debugger, you must first set a breakpoint on an executable statement. Otherwise, the program will run to completion and exit. If the program had a console window, that goes away too. My guess is that this is what is happening.

Try this - select Debug > Start without debugging. What happens?

I ran the example given in

mecej4 — Tue, 01 Mar 2016 09:18:50 GMT

I ran the example given in the IMSL documentation for LINDS. I obtained the correct results and the test code ran with no errors.

As Steve has stated, the normal sequence of events is for the program to run in the debugger to completion, i.e., until a STOP or END statement is executed, if no keyboard input is required. If you did run into difficulties when you ran your code, you will need to provide details of the test case, etc.

When I start without

Pratap__Sangeeta — Tue, 01 Mar 2016 17:58:56 GMT

When I start without debugging, I get a message

*** TERMINAL ERROR 2 FROM LINDS. The required storage cannot be allocated. The workspace is based on N, where N=1258655008.

The matrix I am inverting is 11 by 11.

Steve Lionel (Intel) wrote:

When you run a Fortran program in the Visual Studio debugger, you must first set a breakpoint on an executable statement. Otherwise, the program will run to completion and exit. If the program had a console window, that goes away too. My guess is that this is what is happening.

Try this - select Debug > Start without debugging. What happens?

Sounds as if you passed an

Steven_L_Intel1 — Tue, 01 Mar 2016 18:21:13 GMT

Sounds as if you passed an uninitialized variable. Please show us a small but complete test case.

Here is the program

Pratap__Sangeeta — Tue, 01 Mar 2016 19:22:19 GMT

Here is the program

If I write the whole program in the main program it runs fine. The problem arises when I do it in the subroutine.

Thank you again for your time.

______________

program Test

INCLUDE 'link_fnl_shared.h'
use linds_int

implicit none

integer nf, nt, i, j
parameter (nf=338, nt=11)
real a(nf,nt), apa(nt, nt), apainv(nt, nt), eye(nt, nt)

do i=1,nf
do j=1,nt
a(i,j)=real(i/j)
end do
end do

call getinv(a, nf, nt, apa, apainv)

eye=matmul(apa, apainv)

do i=1, nt
write(*,*) eye(i,:)
pause
end do

  end program Test

subroutine getinv(a, nf, nt, apa, apainv)
external linds

integer nf, nt
real a(nf, nt), ap(nt, nf), apa(nt,nt), apainv(nt, nt)

ap=transpose(a)
apa=matmul(ap, a)

call linds (apa, apainv)


end

You have some rather basic

mecej4 — Wed, 02 Mar 2016 03:14:00 GMT

You have some rather basic errors in your code, and it is to your benefit to read Fortran books and understand why your program did not work.

You call a subroutine (LINDS) that requires an explicit interface, but you did not provide one. The "use linds_int" is not useful in a routine from which LINDS is not called. On the other hand, it must be included in any routine that calls LINDS.

Move the USE LINDS_INT statement from the main program to the subroutine GETINV, remove the EXTERNAL statement, and your program will work.

Thank you for your help. It

Pratap__Sangeeta — Wed, 02 Mar 2016 16:36:05 GMT

Thank you for your help. It does indeed work.

Re. Fortran books, do you have a particular recommendation?

Finally, I find that the computation of the inverse is accurate when everything is declared as double precision (real*8), rather than real. When I declare as real, the inverse multiplied by the matrix was off the identity matrix (diagonal was 1.007 for example, off diagonals were also non zero from the third decimal place ).

In comparison, the Numerical recipes routines ludcmp and lubksb performed well for real variables and I did not have this inaccuracy. .

Is that something one needs to be careful about in all IMSL routines or just this particular one?

Thank you again for your time and patience in answering my rather elementary question.

Best wishes

The Lawrence book is NOT

Steven_L_Intel1 — Wed, 02 Mar 2016 17:37:42 GMT

The Lawrence book is NOT about Fortran - it's how to create Windows API applications in Compaq (and Intel) Visual Fortran.

I mentioned some books in https://software.intel.com/en-us/blogs/2013/12/30/doctor-fortran-in-its-a-modern-fortran-world

John, please mention the

Steven_L_Intel1 — Wed, 02 Mar 2016 19:01:46 GMT

John, please mention the example error in the MKL forum.

Quote:Sangeeta P. wrote:

mecej4 — Thu, 03 Mar 2016 01:20:54 GMT

Sangeeta P. wrote:

I find that the computation of the inverse is accurate when everything is declared as double precision (real*8), rather than real. When I declare as real, the inverse multiplied by the matrix was off the identity matrix (diagonal was 1.007 for example, off diagonals were also non zero from the third decimal place ).

In comparison, the Numerical recipes routines ludcmp and lubksb performed well for real variables and I did not have this inaccuracy. .

Is that something one needs to be careful about in all IMSL routines or just this particular one?

The finding (regarding results with real*8 versus real) is somewhat unexpected. Please present the specific test case in detail. Most of the time, such discrepancies are caused by passing incorrect arguments to routines such as LINDS.

There is always some round

dboggs — Thu, 03 Mar 2016 16:41:53 GMT

There is *always* some round-off error when inverting a matrix using a finite-resolution algorithm. If it is ill-conditioned then the error will show up under a smaller resolution than when it is not. The observation that round-off is serious under REAL*4 than REAL*8 tells me that the reason might be that the subject matrix is ill conditioned, rather than an algorithm error.

Any good text on numerical methods should discuss the conditioning concept.

Dear all

Pratap__Sangeeta — Thu, 03 Mar 2016 16:51:38 GMT

Dear all

Thanks for the comments and suggestions, especially on good reference books.

Here is my program using comparing LINDS and the routines LUDCMP& LUBKSB.

They are not terribly different, but the latter is marginally better. Using REAL and REAL*8 creates a substantial difference. To get the latter, please just change the declarations on the main program and the LUDCMP& LUBKSB subroutines.

Thanks again everyone, for engaging with my problem.

Sangeeta

______________________________________________

program Test

INCLUDE 'link_fnl_shared.h'
use linds_int

implicit none

integer nf, nt, i, j
parameter (nf=338, nt=3)
integer indx(nt)
real a(nf,nt), ap(nt, nf), apa(nt, nt), apainv(nt, nt), eye(nt, nt), b(nf,nt), bp(nt, nf), bpb(nt, nt), bpbinv(nt, nt), eyenew(nt,nt), d, perfecteye(nt, nt), norm_linds, norm_nr

do i=1,nf
do j=1,nt
a(i,j)=real(i/j)
b(i,j)=a(i,j)
end do
end do

ap=transpose(a)
apa=matmul(ap, a)

!Constructing inverse with LINDS
call linds (apa, apainv)
eye=matmul(apa, apainv)

!Constructing inverse with LUDCMP & LUBKSB
bp=transpose(b)
bpb=matmul(bp, b)

!Initializing inverse and creating the accurate identity matrix perfecteye
do i=1, nt
do j=1,nt
bpbinv(i,j)=0.d0
perfecteye(i,j)=0.d0
end do
bpbinv(i,i)=1.d0
perfecteye(i,i)=1.d0
end do

call ludcmp(bpb, nt, nt, indx, d)
do j=1, nt
call lubksb(bpb, nt, nt, indx, bpbinv(1, j))
end do

eyenew=matmul(apa, bpbinv)

!Comparing the distance between perfecteye and the identity matrix got by multiplying the matrix and the inverse
norm_linds=sum(abs(eye-perfecteye))
norm_nr=sum(abs(eyenew-perfecteye))

write(*,*) norm_linds, norm_nr
pause

end program Test

SUBROUTINE ludcmp(a,n,np,indx,d)
INTEGER n,np,indx(n),NMAX
REAL d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0e-20)
INTEGER i,imax,j,k
REAL aamax,dum,sum,vv(NMAX)
d=1.
do 12 i=1,n
aamax=0.
do 11 j=1,n
if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11 continue
if (aamax.eq.0.) pause 'singular matrix in ludcmp'
vv(i)=1./aamax
12 continue
do 19 j=1,n
do 14 i=1,j-1
sum=a(i,j)
do 13 k=1,i-1
sum=sum-a(i,k)*a(k,j)
13 continue
a(i,j)=sum
14 continue
aamax=0.
do 16 i=j,n

sum=a(i,j)
do 15 k=1,j-1
sum=sum-a(i,k)*a(k,j)
15 continue
a(i,j)=sum
dum=vv(i)*abs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
16 continue
if (j.ne.imax)then
do 17 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
17 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.)a(j,j)=TINY
if(j.ne.n)then
dum=1./a(j,j)

do 18 i=j+1,n
a(i,j)=a(i,j)*dum
18 continue
endif
19 continue
return
END
SUBROUTINE lubksb(a,n,np,indx,b)
INTEGER n,np,indx(n)
REAL a(np,np),b(n)
INTEGER i,ii,j,ll
REAL sum
ii=0
do 12 i=1,n
ll=indx(i)
sum=b(ll)
b(ll)=b(i)
if (ii.ne.0)then
do 11 j=ii,i-1
sum=sum-a(i,j)*b(j)
11 continue
else if (sum.ne.0.) then
ii=i
endif
b(i)=sum
12 continue
do 14 i=n,1,-1
sum=b(i)
do 13 j=i+1,n
sum=sum-a(i,j)*b(j)
13 continue
b(i)=sum/a(i,i)
14 continue
return
END

I have this comment regarding

mecej4 — Fri, 04 Mar 2016 03:40:09 GMT

I have this comment regarding the first part of the program code in #20.

You compute X = A^TA for a matrix A of size 338 X 3 and then find the inverse of X. For the particular matrix that you chose, i.e., a_ij = floor(i/j), the condition number of A is about 930, and the condition number of X is the square of that, i.e., about 9 X 10⁵. If you calculate X in single precision, that means that only one or two significant digits may be correct.

Two further points to note: in general, it is bad to form the normal matrix X and then work with it. If your ultimate objective is to solve overdetermined linear equations A.x = b, for example, rather than forming X, use the well-known QR decomposition of A. Likewise, it is bad practice to form the inverse matrix explictly if your objective is to solve linear equations. When I say "bad", I mean (i) inefficient and (ii) inaccurate. Any book on numerical analysis and Wikipedia articles can provide you with details on why the "bad" happens.

Fortran Books:

JohnNichols — Wed, 02 Mar 2016 17:30:45 GMT

Fortran Books:

Lawrences - Compaq Visual Fortran is not bad - a bit dated

The Metcalf and Reid books are all good - a bit to academic at times.

If you want to learn to program well - try Abelson and Sussman or Winston and Horne - all LISP - but if you LISP well - you are a long way forward.

The old MICROSOFT Fortran manuals are wonderful if you can get them.

JohnNichols — Wed, 02 Mar 2016 17:37:11 GMT

In comparison, the Numerical recipes routines ludcmp and lubksb performed well for real variables and I did not have this inaccuracy. .

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

I have been playng with three inverters, one from Harrison in 1973, one from Conte and deBoore which is mid 70's and I just got through getting PARDISO to run with these problems.

All gave the inversions to an accuracy of at least six figures and agreed.

!===============================================================================
! Copyright 2004-2016 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

!===============================================================================
! Copyright 2004-2016 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

*   Content : MKL PARDISO Fortran example
*
********************************************************************************
C----------------------------------------------------------------------
C Example program to show the use of the "PARDISO" routine
C for symmetric linear systems
C---------------------------------------------------------------------
C This program can be downloaded from the following site:
C www.pardiso-project.org
C
C (C) Olaf Schenk, Department of Computer Science,
C University of Basel, Switzerland.
C Email: olaf.schenk@unibas.ch
C
C---------------------------------------------------------------------

      MODULE Solver
      INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
      REAL (KIND=dp) :: g = 9.806, pi = 3.14159265D0

      contains

      subroutine linsolve(n, error, nnz, a, ja, ia, b)
      implicit none
      include 'mkl_pardiso.fi'

      INTEGER :: n, nnz
      REAL (KIND=dp) :: x(n) ! Pardiso needs X(N) even if sol is returned in b
      
      REAL (KIND=dp) :: b(n)
C.. Internal solver memory pointer for 64-bit architectures
C.. INTEGER*8 pt(64)
C.. Internal solver memory pointer for 32-bit architectures
C.. INTEGER*4 pt(64)
C.. This is OK in both cases
      TYPE(MKL_PARDISO_HANDLE) pt(64)
C.. All other variables
      INTEGER maxfct, mnum, mtype, phase, nrhs, msglvl
      INTEGER iparm(64)
      INTEGER ia(nnz)
      INTEGER ja(nnz)
      INTEGER i, idum(1)
      REAL (KIND=dp) ::  ddum(1)
      integer error
      REAL (KIND=dp) a(nnz)
C.. Fill all arrays containing matrix data.
      DATA nrhs /1/, maxfct /1/, mnum /1/
C..
C.. Setup PARDISO control parameter
C..
      DO i = 1, 64
      iparm(i) = 0
      END DO
      iparm(1) = 1 ! no solver default
      iparm(2) = 2 ! fill-in reordering from METIS
      iparm(3) = 1 ! numbers of processors
      iparm(4) = 0 ! no iterative-direct algorithm
      iparm(5) = 0 ! no user fill-in reducing permutation
      iparm(6) = 0 ! =0 solution on the first n components of x
      iparm(7) = 0 ! not in use
      iparm(8) = 9 ! numbers of iterative refinement steps
      iparm(9) = 0 ! not in use
      iparm(10) = 13 ! perturb the pivot elements with 1E-13
      iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
      iparm(12) = 0 ! not in use
      iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
      iparm(14) = 0 ! Output: number of perturbed pivots
      iparm(15) = 0 ! not in use
      iparm(16) = 0 ! not in use
      iparm(17) = 0 ! not in use
      iparm(18) = -1 ! Output: number of nonzeros in the factor LU
      iparm(19) = -1 ! Output: Mflops for LU factorization
      iparm(20) = 0 ! Output: Numbers of CG Iterations
      error = 0 ! initialize error flag
      msglvl = 1 ! print statistical information
      mtype = 11 ! real unsymmetric
C.. Initialize the internal solver memory pointer. This is only
C necessary for the FIRST call of the PARDISO solver.
      DO i = 1, 64
          pt(i)%DUMMY = 0
      END DO
C.. Reordering and Symbolic Factorization, This step also allocates
C   all memory that is necessary for the factorization
      phase = 11 ! only reordering and symbolic factorization
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     &idum, nrhs, iparm, msglvl, ddum, ddum, error)
      WRITE(*,*) 'Reordering completed ... '
      IF (error .NE. 0) THEN
          WRITE(*,*) 'The following ERROR was detected: ', error
          STOP 1
      END IF
      WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
      WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
C.. Factorization.
      phase = 22 ! only factorization
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     &idum, nrhs, iparm, msglvl, ddum, ddum, error)
      WRITE(*,*) 'Factorization completed ... '
      IF (error .NE. 0) THEN
          WRITE(*,*) 'The following ERROR was detected: ', error
          STOP 1
      END IF
C.. Back substitution and iterative refinement
      iparm(8) = 2 ! max numbers of iterative refinement steps
      phase = 33 ! only factorization
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     &idum, nrhs, iparm, msglvl, b, x, error)
      WRITE(*,*) 'Solve completed ... '
      WRITE(*,*) 'The solution of the system is '
      DO i = 1, n
          WRITE(*,*) ' x(',i,') = ', x(i)
      END DO
C.. Termination and release of memory
      phase = -1 ! release internal memory
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum,
     &idum, idum, nrhs, iparm, msglvl, ddum, ddum, error)
      Return
      End subroutine linsolve

      end module

The driver routine for

JohnNichols — Wed, 02 Mar 2016 17:48:33 GMT

The driver routine for PARDISO to take a square matrix and put it into the correct form is shown - it is not pretty - merely allows me to test LINSOLVE before I use it for real.

ludcmp and lubksb -- PARDISO is a lot faster -- a.txt is the sample from the INTEL examples.

Lawrence is a good book --

!------------------------------------------------------------------------------------------------------------------------
!     
!     Test Problem
!
!------------------------------------------------------------------------------------------------------------------------
      
!      ja(1) = 1
!      ja(2) = 2
!      ja(3) = 4
!      ja(4) = 1
!      ja(5) = 2
!      ja(6) = 3
!      ja(7) = 4
!      ja(8) = 5
!      ja(9) = 1
!      ja(10) = 3
!      ja(11) = 4
!      ja(12) = 2
!      ja(13) = 5

!      ia(1) = 1
!      ia(2) = 4
!      ia(3) = 6
!      ia(4) = 9
!      ia(5) = 12
!      ia(6) = 14

!      a(1) = 1.d0
!      a(2) = -1.d0
!      a(3) = -3.d0
!      a(4) = -2.d0
!      a(5) = 5.d0
!      a(6) = 4.d0
!      a(7) = 6.d0
!      a(8) = 4.d0
!      a(9) = -4.d0
!      a(10)= 2.d0
!      a(11) = -7.d0
!      a(12) = 8.d0
!      a(13) = -5.d0
      
!DATA a
! 1  / 1.d0,-1.d0,      -3.d0,
! 2   -2.d0, 5.d0,
!  3                4.d0, 6.d0, 4.d0,
!   4   -4.d0,       2.d0, 7.d0,
!   5          8.d0,            -5.d0/
      
      
      PROGRAM pardiso_unsym

      use Solver
      IMPLICIT NONE

      include 'mkl_pardiso.fi'


      REAL (KIND=dp), ALLOCATABLE :: a(:), b(:), c(:,:), atemp(:)
      integer, allocatable :: ja(:), ia(:), iatemp(:), jatemp(:)
      INTEGER :: nnodes = 5
      integer error, iaN
      integer :: nnz = 200, istat, i, j , count
      integer :: size = 0
      integer :: size1 = 0
      integer :: size2 = 0
      integer :: flag = 0
      REAL (KIND=dp) :: delta = 0.000000000000001


      open(1, file="a.txt", STATUS = 'old')

      read(1,*)count
      
      nnodes = count      
      iaN = count+1
      
      ALLOCATE (a(nnz), ja(nnz), ia(iaN), jatemp(nnz), iatemp(iaN), b(count), c(count,count), atemp(nnz), STAT=istat)
      IF (istat.NE.0) THEN
            WRITE (*, *) '*** Could not allocate some arrays in LINSOLVE'
            STOP
      END IF

      b = 1.0D0
      a = 0.0d0
      c = 0.0d0
      iatemp = 0
      jatemp = 0
      atemp = 0
      ia = 0
      ja = 0
      do 200 i = 1, count
            flag = 0
            do 300 j = 1, count
                  read(1,*)c(i,j)
                  if(abs(c(i,j)) .gt. delta) then
                        write(*,*)c(i,j)
                        size = size + 1
                        jatemp(size) = j
                        iatemp(size1 + 1) = size+1
                        atemp(size) = c(i,j)
                        if(flag .eq. 0) then
                              size1 = size1 + 1
                              iatemp(size1) = size

                              write(*,*)size1, iatemp(size1)
                              flag = 1
                        endif
                  endif
300   end do

      !    write(*,400)(c(i,j),j=1,count)
400   format(5(1x,F10.3))
200   end do

      do 500 i = 1, count
            read(1,*)b(i)
500   end  do
      write(*,100)count, size
100   Format(' Matrix Count :: ',i5, ' Size :: ',i5)


      call LinSolve(nnodes, nnz, error, atemp, jatemp, iatemp, b)

      end

for symmetric linear systems

JohnNichols — Wed, 02 Mar 2016 17:51:35 GMT

for symmetric linear systems

Steve:

There is a small error in the example for the non-symmetric matrix in the INTEL examples - it says symmetric in the FORTRAN text file.

John