Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Cannot invert a matrix with IMSL routine LINDS

Pratap__Sangeeta
Beginner
1,094 Views

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
TimP
Honored Contributor III
992 Views

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
Pratap__Sangeeta
Beginner
992 Views

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

Thanks!

0 Kudos
Steven_L_Intel1
Employee
992 Views

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
mecej4
Honored Contributor III
992 Views

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
Pratap__Sangeeta
Beginner
992 Views

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
Steven_L_Intel1
Employee
992 Views

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

0 Kudos
Pratap__Sangeeta
Beginner
992 Views

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
mecej4
Honored Contributor III
992 Views

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
Pratap__Sangeeta
Beginner
992 Views

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
JohnNichols
Valued Contributor III
990 Views

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
JohnNichols
Valued Contributor III
990 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
Steven_L_Intel1
Employee
992 Views

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

0 Kudos
JohnNichols
Valued Contributor III
990 Views

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
JohnNichols
Valued Contributor III
990 Views

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
Steven_L_Intel1
Employee
992 Views

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

0 Kudos
mecej4
Honored Contributor III
992 Views

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
JohnNichols
Valued Contributor III
990 Views

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
dboggs
New Contributor I
992 Views

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
Pratap__Sangeeta
Beginner
992 Views

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
JohnNichols
Valued Contributor III
805 Views

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

your inversion matrix is so dense and not banded that it is not funny - this is your a matrix run through PARDISO - adjust your problem to something that is realistic and maybe it will work as you want.

   1.00000000000000
           1           1
   2.00000000000000
           2           2
   1.00000000000000
   3.00000000000000
           3           4
   1.00000000000000
   1.00000000000000
   4.00000000000000
           4           7
   2.00000000000000
   1.00000000000000
   1.00000000000000
   5.00000000000000
           5          11
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   6.00000000000000
           6          16
   3.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   7.00000000000000
           7          22
   3.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   8.00000000000000
           8          29
   4.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   9.00000000000000
           9          37
   4.00000000000000
   3.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   10.0000000000000
          10          46
   5.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   11.0000000000000
          11          56
   5.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   12.0000000000000
          12          67
   6.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   13.0000000000000
          13          79
   6.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   14.0000000000000
          14          92
   7.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   15.0000000000000
          15         106
   7.00000000000000
   5.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   16.0000000000000
          16         121
   8.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   17.0000000000000
          17         137
   8.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   18.0000000000000
          18         154
   9.00000000000000
   6.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   19.0000000000000
          19         172
   9.00000000000000
   6.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   20.0000000000000
          20         191
   10.0000000000000
   6.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   21.0000000000000
          21         211
   10.0000000000000
   7.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   22.0000000000000
          22         232
   11.0000000000000
   7.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   23.0000000000000
          23         254
   11.0000000000000
   7.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   24.0000000000000
          24         277
   12.0000000000000
   8.00000000000000
   6.00000000000000
   4.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   25.0000000000000
          25         301
   12.0000000000000
   8.00000000000000
   6.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   26.0000000000000
          26         326
   13.0000000000000
   8.00000000000000
   6.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   27.0000000000000
          27         352
   13.0000000000000
   9.00000000000000
   6.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   28.0000000000000
          28         379
   14.0000000000000
   9.00000000000000
   7.00000000000000
   5.00000000000000
   4.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   29.0000000000000
          29         407
   14.0000000000000
   9.00000000000000
   7.00000000000000
   5.00000000000000
   4.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   30.0000000000000
          30         436
   15.0000000000000
   10.0000000000000
   7.00000000000000
   6.00000000000000
   5.00000000000000
   4.00000000000000
   3.00000000000000
   3.00000000000000
   3.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   2.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
   1.00000000000000
 Matrix Count ::    30 Size ::   465

=== PARDISO: solving a real nonsymmetric system ===
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Scaling is turned ON
Matching is turned ON


Summary: ( reordering phase )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000008 s
Time spent in reordering of the initial matrix (reorder)         : 0.000149 s
Time spent in symbolic factorization (symbfct)                   : 0.000755 s
Time spent in data preparations for factorization (parlist)      : 0.000001 s
Time spent in allocation of internal data structures (malloc)    : 0.001801 s
Time spent in additional calculations                            : 0.000024 s
Total time spent                                                 : 0.002738 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           30
             number of non-zeros in A:      465
             number of non-zeros in A (%): 51.666667

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    1
             size of largest supernode:               30
             number of non-zeros in L:                900
             number of non-zeros in U:                1
             number of non-zeros in L+U:              901
 Reordering completed ...
 Number of nonzeros in factors =          901
 Number of factorization MFLOPS =            0
=== PARDISO is running in In-Core mode, because iparam(60)=0 ===

Percentage of computed non-zeros for LL^T factorization
 100 %


=== PARDISO: solving a real nonsymmetric system ===
Single-level factorization algorithm is turned ON


Summary: ( factorization phase )
================

Times:
======
Time spent in copying matrix to internal data structure (A to LU): 0.000000 s
Time spent in factorization step (numfct)                        : 0.001749 s
Time spent in allocation of internal data structures (malloc)    : 0.001294 s
Time spent in additional calculations                            : 0.000004 s
Total time spent                                                 : 0.003048 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           30
             number of non-zeros in A:      465
             number of non-zeros in A (%): 51.666667

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    1
             size of largest supernode:               30
             number of non-zeros in L:                900
             number of non-zeros in U:                1
             number of non-zeros in L+U:              901
             gflop   for the numerical factorization: 0.000018

             gflop/s for the numerical factorization: 0.010048

 Factorization completed ...

=== PARDISO: solving a real nonsymmetric system ===


Summary: ( solution phase )
================

Times:
======
Time spent in direct solver at solve step (solve)                : 0.000239 s
Time spent in additional calculations                            : 0.000730 s
Total time spent                                                 : 0.000969 s

Statistics:
===========
Parallel Direct Factorization is running on 4 OpenMP

< Linear system Ax = b >
             number of equations:           30
             number of non-zeros in A:      465
             number of non-zeros in A (%): 51.666667

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with state of the art partitioning metis>
             number of supernodes:                    1
             size of largest supernode:               30
             number of non-zeros in L:                900
             number of non-zeros in U:                1
             number of non-zeros in L+U:              901
             gflop   for the numerical factorization: 0.000018

             gflop/s for the numerical factorization: 0.010048

 Solve completed ...
 The solution of the system is
  x(           1 ) =   0.000000000000000E+000
  x(           2 ) =   -21.0000000000000
  x(           3 ) =    21.0000000000000
  x(           4 ) =    21.0000000000000
  x(           5 ) =   1.878240250526247E-032
  x(           6 ) =   7.771561172376110E-017
  x(           7 ) =  -1.436853791652556E-031
  x(           8 ) =   -21.0000000000000
  x(           9 ) =  -7.888609052210118E-031
  x(          10 ) =    21.0000000000000
  x(          11 ) =   8.545993139894292E-031
  x(          12 ) =   -21.0000000000000
  x(          13 ) =   2.253888300631437E-032
  x(          14 ) =    21.0000000000000
  x(          15 ) =   -21.0000000000000
  x(          16 ) =    21.0000000000000
  x(          17 ) =   0.000000000000000E+000
  x(          18 ) =   7.027711745877242E-015
  x(          19 ) =   0.000000000000000E+000
  x(          20 ) =   -21.0000000000000
  x(          21 ) =   -21.0000000000000
  x(          22 ) =    21.0000000000000
  x(          23 ) =   0.000000000000000E+000
  x(          24 ) =    21.0000000000000
  x(          25 ) =   1.972152263052530E-031
  x(          26 ) =    21.0000000000000
  x(          27 ) =   -21.0000000000000
  x(          28 ) =   -21.0000000000000
  x(          29 ) =   0.000000000000000E+000
  x(          30 ) =  -5.995204332975852E-015
Press any key to continue . . .

0 Kudos
Reply