稀疏矩阵高性能求解器PARDISO的使用
INCLUDE 'mkl_pardiso.f90'PROGRAM pardiso_sym_f90USE mkl_pardisoIMPLICIT NONEINTEGER, PARAMETER :: dp = KIND(1.0D0)!.. Internal solver memory pointerTYPE(MKL_PARDISO_HANDLE), ALLOCATABLE :: pt(:)!.. All other variablesINTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl, nnzINTEGER error1INTEGER, ALLOCATABLE :: iparm( : )INTEGER, ALLOCATABLE :: ia( : )INTEGER, ALLOCATABLE :: ja( : )REAL(KIND=DP), ALLOCATABLE :: a( : )REAL(KIND=DP), ALLOCATABLE :: b( : )REAL(KIND=DP), ALLOCATABLE :: x( : )INTEGER i, idum(1)REAL(KIND=DP) ddum(1)!.. Fill all arrays containing matrix data.n = 8nnz = 18nrhs = 1maxfct = 1mnum = 1ALLOCATE(ia(n + 1))ia = (/ 1, 5, 8, 10, 12, 15, 17, 18, 19 /)ALLOCATE(ja(nnz))ja = (/ 1, 3, 6, 7, &2, 3, 5, &3, 8, &4, 7, &5, 6, 7, &6, 8, &7, &8 /)ALLOCATE(a(nnz))a = (/ 7.d0, 1.d0, 2.d0, 7.d0, &-4.d0, 8.d0, 2.d0, &1.d0, 5.d0, &7.d0, 9.d0, &5.d0, 1.d0, 5.d0, &-1.d0, 5.d0, &11.d0, &5.d0 /)ALLOCATE(b(n))ALLOCATE(x(n))!..!.. Set up PARDISO control parameter!..ALLOCATE(iparm(64))DO i = 1, 64iparm(i) = 0END DOiparm(1) = 1 ! no solver defaultiparm(2) = 2 ! fill-in reordering from METISiparm(4) = 0 ! no iterative-direct algorithmiparm(5) = 0 ! no user fill-in reducing permutationiparm(6) = 0 ! =0 solution on the first n components of xiparm(8) = 2 ! numbers of iterative refinement stepsiparm(10) = 13 ! perturb the pivot elements with 1E-13iparm(11) = 1 ! use nonsymmetric permutation and scaling MPSiparm(13) = 0iparm(14) = 0 ! Output: number of perturbed pivotsiparm(18) = -1 ! Output: number of nonzeros in the factor LUiparm(19) = -1 ! Output: Mflops for LU factorizationiparm(20) = 0 ! Output: Numbers of CG Iterationserror = 0 ! initialize error flagmsglvl = 1 ! print statistical informationmtype = -2 ! symmetric, indefinite!.. Initialize the internal solver memory pointer. This is only! necessary for the FIRST call of the PARDISO solver.ALLOCATE (pt(64))DO i = 1, 64pt(i)%DUMMY = 0END DO!.. Reordering and Symbolic Factorization,! This step also allocates! all memory that is necessary for the factorizationphase = 11 ! only reordering and symbolic factorizationCALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &idum, nrhs, iparm, msglvl, ddum, ddum, error)WRITE(*,*) 'Reordering completed ... 'IF (error /= 0) THENWRITE(*,*) 'The following ERROR was detected: ', errorGOTO 1000END IFWRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)!.. Factorization.phase = 22 ! only factorizationCALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &idum, nrhs, iparm, msglvl, ddum, ddum, error)WRITE(*,*) 'Factorization completed ... 'IF (error /= 0) THENWRITE(*,*) 'The following ERROR was detected: ', errorGOTO 1000ENDIF!.. Back substitution and iterative refinementiparm(8) = 2 ! max numbers of iterative refinement stepsphase = 33 ! only solvingDO i = 1, nb(i) = 1.d0END DOCALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &idum, nrhs, iparm, msglvl, b, x, error)WRITE(*,*) 'Solve completed ... 'IF (error /= 0) THENWRITE(*,*) 'The following ERROR was detected: ', errorGOTO 1000ENDIFWRITE(*,*) 'The solution of the system is 'DO i = 1, nWRITE(*,*) ' x(',i,') = ', x(i)END DO1000 CONTINUE!.. Termination and release of memoryphase = -1 ! release internal memoryCALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &idum, nrhs, iparm, msglvl, ddum, ddum, error1)IF (ALLOCATED(ia)) DEALLOCATE(ia)IF (ALLOCATED(ja)) DEALLOCATE(ja)IF (ALLOCATED(a)) DEALLOCATE(a)IF (ALLOCATED(b)) DEALLOCATE(b)IF (ALLOCATED(x)) DEALLOCATE(x)IF (ALLOCATED(iparm)) DEALLOCATE(iparm)IF (error1 /= 0) THENWRITE(*,*) 'The following ERROR on release stage was detected: ', error1STOP 1ENDIFIF (error /= 0) STOP 1END PROGRAM pardiso_sym_f90
// Define the format to printf MKL_INT valuesint main(void){/* Matrix data. */MKL_INT n = 8;MKL_INT ia[9] = { 0, 4, 7, 9, 11, 14, 16, 17, 18 };MKL_INT ja[18] ={ 0, 2, 5, 6,1, 2, 4,2, 7,3, 6,4, 5, 6,5, 7,6,7};double a[18] ={ 7.0, 1.0, 2.0, 7.0,-4.0, 8.0, 2.0,1.0, 5.0,7.0, 9.0,5.0, 1.0, 5.0,-1.0, 5.0,11.0,5.0};MKL_INT mtype = -2; /* Real symmetric matrix *//* RHS and solution vectors. */double b[8], x[8];MKL_INT nrhs = 1; /* Number of right hand sides. *//* Internal solver memory pointer pt, *//* 32-bit: int pt[64]; 64-bit: long int pt[64] *//* or void *pt[64] should be OK on both architectures */void* pt[64];/* Pardiso control parameters. */MKL_INT iparm[64];MKL_INT maxfct, mnum, phase, error, msglvl;/* Auxiliary variables. */MKL_INT i;double ddum; /* Double dummy */MKL_INT idum; /* Integer dummy. *//* -------------------------------------*//* .. Setup Pardiso control parameters. *//* -------------------------------------*/for (i = 0; i < 64; i++){iparm[i] = 0;}iparm[0] = 1; /* No solver default */iparm[1] = 2; /* Fill-in reordering from METIS */iparm[3] = 0; /* No iterative-direct algorithm */iparm[4] = 0; /* No user fill-in reducing permutation */iparm[5] = 0; /* Write solution into x */iparm[7] = 2; /* Max numbers of iterative refinement steps */iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */iparm[12] = 0;iparm[13] = 0; /* Output: Number of perturbed pivots */iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */iparm[18] = -1; /* Output: Mflops for LU factorization */iparm[19] = 0; /* Output: Numbers of CG Iterations */iparm[34] = 1; /* PARDISO use C-style indexing for ia and ja arrays */maxfct = 1; /* Maximum number of numerical factorizations. */mnum = 1; /* Which factorization to use. */msglvl = 1; /* Print statistical information in file */error = 0; /* Initialize error flag *//* ----------------------------------------------------------------*//* .. Initialize the internal solver memory pointer. This is only *//* necessary for the FIRST call of the PARDISO solver. *//* ----------------------------------------------------------------*/for (i = 0; i < 64; i++){pt[i] = 0;}/* --------------------------------------------------------------------*//* .. Reordering and Symbolic Factorization. This step also allocates *//* all memory that is necessary for the factorization. *//* --------------------------------------------------------------------*/phase = 11;PARDISO(pt, &maxfct, &mnum, &mtype, &phase,&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);if (error != 0){printf("\nERROR during symbolic factorization: " IFORMAT, error);exit(1);}printf("\nReordering completed ... ");printf("\nNumber of nonzeros in factors = " IFORMAT, iparm[17]);printf("\nNumber of factorization MFLOPS = " IFORMAT, iparm[18]);/* ----------------------------*//* .. Numerical factorization. *//* ----------------------------*/phase = 22;PARDISO(pt, &maxfct, &mnum, &mtype, &phase,&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);if (error != 0){printf("\nERROR during numerical factorization: " IFORMAT, error);exit(2);}printf("\nFactorization completed ... ");/* -----------------------------------------------*//* .. Back substitution and iterative refinement. *//* -----------------------------------------------*/phase = 33;iparm[7] = 2; /* Max numbers of iterative refinement steps. *//* Set right hand side to one. */for (i = 0; i < n; i++){b[i] = 1;}PARDISO(pt, &maxfct, &mnum, &mtype, &phase,&n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);if (error != 0){printf("\nERROR during solution: " IFORMAT, error);exit(3);}printf("\nSolve completed ... ");printf("\nThe solution of the system is: ");for (i = 0; i < n; i++){printf("\n x [" IFORMAT "] = % f", i, x[i]);}printf("\n");/* --------------------------------------*//* .. Termination and release of memory. *//* --------------------------------------*/phase = -1; /* Release internal memory. */PARDISO(pt, &maxfct, &mnum, &mtype, &phase,&n, &ddum, ia, ja, &idum, &nrhs,iparm, &msglvl, &ddum, &ddum, &error);return 0;}
登录后免费查看全文
附件
免费链接.txt
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2023-07-17
最近编辑：1年前
硕士
|
cae工程师
签名征集中
相关推荐
热门文章
