Community
cancel
Showing results for 
Search instead for 
Did you mean: 
JohnNichols
Valued Contributor I
48 Views

Error in the comment on non-symmetric PARDISO - should say non-

!===============================================================================
! 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

This is called symmetirc in the comments - but it comes from the non-symmetric file example.

This is a simple example to take a square matrix and turn it into a PARDISO style matricx with zeros removed.  -  The error can be seen in the code comments. A.txt is the intel sample data in a square format

!------------------------------------------------------------------------------------------------------------------------
!     
!     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.mat", 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
1 Reply
Gennady_F_Intel
Moderator
48 Views

ok, thanks for the case

Reply