Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
2 Views

Cannot invert a matrix with IMSL routine LINDS

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

 

0 Kudos
22 Replies
Highlighted
Black Belt
2 Views

This would likely happen if

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.

0 Kudos
Highlighted
2 Views

No, I am using Intel Visual

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

Thanks!

0 Kudos
Highlighted
2 Views

When you run a Fortran

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?

Retired 12/31/2016
0 Kudos
Highlighted
Black Belt
2 Views

I ran the example given in

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. 

0 Kudos
Highlighted
2 Views

When I start without

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?

0 Kudos
Highlighted
2 Views

Sounds as if you passed an

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

Retired 12/31/2016
0 Kudos
Highlighted
2 Views

Here is the program

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

0 Kudos
Highlighted
Black Belt
2 Views

You have some rather basic

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.

 

0 Kudos
Highlighted
2 Views

Thank you for your help. It

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

 

0 Kudos
Highlighted
New Contributor I
2 Views

Fortran Books:

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. 

 

0 Kudos
Highlighted
New Contributor I
2 Views

 

 

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

 

0 Kudos
Highlighted
2 Views

The Lawrence book is NOT

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

Retired 12/31/2016
0 Kudos
Highlighted
New Contributor I
2 Views

The driver routine for

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

 

0 Kudos
Highlighted
New Contributor I
2 Views

for symmetric linear systems

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

0 Kudos
Highlighted
2 Views

John, please mention the

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

Retired 12/31/2016
0 Kudos
Highlighted
Black Belt
2 Views

Quote:Sangeeta P. wrote:

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.

0 Kudos
Highlighted
New Contributor I
2 Views

Quote:mecej4 wrote:

mecej4 wrote:

Quote:

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.

Given the pure mathematics behind the different inversion routines and the long history of their use from the 1950's I would not believe the original statement -- I agree with mecej4, one needs to present a solid working example to make such an assertion and then the error can be located.  I doubt it is an error in the supplied rotuines.

0 Kudos
Highlighted
New Contributor I
2 Views

There is *always* some round

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. 

0 Kudos
Highlighted
2 Views

Dear all

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

0 Kudos