- 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