- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
!===============================================================================
! 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
Link Copied
1 Reply
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
ok, thanks for the case
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page