/* -------------------------------------------------------------------- */ /* Example program to show the use of the "PARDISO" routine */ /* on for unsymmetric linear systems */ /* -------------------------------------------------------------------- */ /* This program can be downloaded from the following site: */ /* http://www.panua.ch/products/pardiso/ */ /* */ /* (C) Olaf Schenk, Panua Technologies, Switzerland. */ /* Email: info@panua.ch */ /* -------------------------------------------------------------------- */ #include #include #include /* PARDISO prototype. */ void pardisoinit (void *, int *, int *, int *, double *, int *); void pardiso (void *, int *, int *, int *, int *, int *, double *, int *, int *, int *, int *, int *, int *, double *, double *, int *, double *); void pardiso_chkmatrix (int *, int *, double *, int *, int *, int *); void pardiso_chkvec (int *, int *, double *, int *); void pardiso_printstats (int *, int *, double *, int *, int *, int *, double *, int *); int main( void ) { /* Matrix data. */ int n = 8; int ia[ 9] = { 0, 4, 7, 9, 11, 12, 15, 17, 20 }; int ja[20] = { 0, 2, 5, 6, 1, 2, 4, 2, 7, 3, 6, 1, 2, 5, 7, 1, 6, 2, 6, 7 }; double a[20] = { 7.0, 1.0, 2.0, 7.0, -4.0, 8.0, 2.0, 1.0, 5.0, 7.0, 9.0, -4.0, 7.0, 3.0, 8.0, 1.0, 11.0, -3.0, 2.0, 5.0 }; int nnz = ia[n]; int mtype = 11; /* Real unsymmetric matrix */ /* RHS and solution vectors. */ double b[8], x[8], diag[8]; 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. */ int iparm[64]; double dparm[64]; int solver; int maxfct, mnum, phase, error, msglvl; /* Number of processors. */ int num_procs; /* Auxiliary variables. */ char *var; int i, k; double ddum; /* Double dummy */ int idum; /* Integer dummy. */ /* -------------------------------------------------------------------- */ /* .. Setup Pardiso control parameters and initialize the solvers */ /* internal adress pointers. This is only necessary for the FIRST */ /* call of the PARDISO solver. */ /* ---------------------------------------------------------------------*/ error = 0; solver = 0; /* use sparse direct solver */ pardisoinit (pt, &mtype, &solver, iparm, dparm, &error); if (error != 0) { if (error == -10 ) printf("No license file found \n"); if (error == -11 ) printf("License is expired \n"); if (error == -12 ) printf("Wrong username or hostname \n"); return 1; } else printf("[PARDISO]: License check was successful ... \n"); /* Numbers of processors, value of OMP_NUM_THREADS */ var = getenv("OMP_NUM_THREADS"); if(var != NULL) sscanf( var, "%d", &num_procs ); else { printf("Set environment OMP_NUM_THREADS to 1"); exit(1); } iparm[2] = num_procs; iparm[10] = 1; /* no scaling */ iparm[12] = 1; /* no matching */ iparm[26] = 1; /* recaling */ maxfct = 1; /* Maximum number of numerical factorizations. */ mnum = 1; /* Which factorization to use. */ msglvl = 1; /* Print statistical information */ error = 0; /* Initialize error flag */ /* -------------------------------------------------------------------- */ /* .. Convert matrix from 0-based C-notation to Fortran 1-based */ /* notation. */ /* -------------------------------------------------------------------- */ for (i = 0; i < n+1; i++) { ia[i] += 1; } for (i = 0; i < nnz; i++) { ja[i] += 1; } /* Set right hand side to i. */ for (i = 0; i < n; i++) { b[i] = i; } /* -------------------------------------------------------------------- */ /* .. pardiso_chk_matrix(...) */ /* Checks the consistency of the given matrix. */ /* Use this functionality only for debugging purposes */ /* -------------------------------------------------------------------- */ pardiso_chkmatrix (&mtype, &n, a, ia, ja, &error); if (error != 0) { printf("\nERROR in consistency of matrix: %d", error); exit(1); } /* -------------------------------------------------------------------- */ /* .. pardiso_chkvec(...) */ /* Checks the given vectors for infinite and NaN values */ /* Input parameters (see PARDISO user manual for a description): */ /* Use this functionality only for debugging purposes */ /* -------------------------------------------------------------------- */ pardiso_chkvec (&n, &nrhs, b, &error); if (error != 0) { printf("\nERROR in right hand side: %d", error); exit(1); } /* -------------------------------------------------------------------- */ /* .. pardiso_printstats(...) */ /* prints information on the matrix to STDOUT. */ /* Use this functionality only for debugging purposes */ /* -------------------------------------------------------------------- */ pardiso_printstats (&mtype, &n, a, ia, ja, &nrhs, b, &error); if (error != 0) { printf("\nERROR right hand side: %d", error); exit(1); } /* -------------------------------------------------------------------- */ /* .. 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, dparm); if (error != 0) { printf("\nERROR during symbolic factorization: %d", error); exit(1); } printf("\nReordering completed ... "); printf("\nNumber of nonzeros in factors = %d", iparm[17]); printf("\nNumber of factorization GFLOPS = %d", iparm[18]); /* -------------------------------------------------------------------- */ /* .. Numerical factorization. */ /* -------------------------------------------------------------------- */ phase = 22; pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error, dparm); if (error != 0) { printf("\nERROR during numerical factorization: %d", error); exit(2); } printf("\nFactorization completed ...\n "); /* -------------------------------------------------------------------- */ /* .. Numerical factorization. */ /* -------------------------------------------------------------------- */ phase = 22; a[0] = 100*a[0]; pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error, dparm); if (error != 0) { printf("\nERROR during numerical factorization: %d", error); exit(2); } printf("\nFactorization completed ...\n "); /* -------------------------------------------------------------------- */ /* .. Back substitution and iterative refinement. */ /* -------------------------------------------------------------------- */ phase = 33; iparm[7] = 1; /* Max numbers of iterative refinement steps. */ pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm); if (error != 0) { printf("\nERROR during solution: %d", error); exit(3); } printf("\nSolve completed ... "); printf("\nThe solution of the system is: "); for (i = 0; i < n; i++) { printf("\n x [%d] = % f", i, x[i] ); } printf ("\n"); /* -------------------------------------------------------------------- */ /* .. Back substitution with tranposed matrix A^t x=b */ /* -------------------------------------------------------------------- */ phase = 33; iparm[7] = 1; /* Max numbers of iterative refinement steps. */ iparm[11] = 1; /* Solving with transpose matrix. */ pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm); if (error != 0) { printf("\nERROR during solution: %d", error); exit(3); } printf("\nSolve completed ... "); printf("\nThe solution of the system is: "); for (i = 0; i < n; i++) { printf("\n x [%d] = % f", i, x[i] ); } printf ("\n"); /* -------------------------------------------------------------------- */ /* ... compute diagonal elements of the inverse. */ /* -------------------------------------------------------------------- */ phase = 33; iparm[11] = 0; /* Solving with nontranspose matrix. */ /* solve for n right hand sides */ for (k = 0; k < n; k++) { for (i = 0; i < n; i++) { b[i] = 0; } /* Set k-th right hand side to one. */ b[k] = 1; pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm); if (error != 0) { printf("\nERROR during solution: %d", error); exit(3); } /* save diagonal element */ diag[k] = x[k]; } /* -------------------------------------------------------------------- */ /* ... Inverse factorization. */ /* -------------------------------------------------------------------- */ if (solver == 0) { printf("\nCompute Diagonal Elements of the inverse of A ... \n"); phase = -22; iparm[35] = 0; /* overwrite internal factor L */ pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error, dparm); /* print diagonal elements */ for (i = 0; i < n; i++) { for (int j = ia[i]-1; j < ia[i+1]-1; j++) { int k = ja[j]-1; if (k == i) { printf ("Diagonal element %d of A^{-1} = %32.24e = %32.24e \n", i, a[j], diag[i]); break; } } } } /* -------------------------------------------------------------------- */ /* .. Convert matrix back to 0-based C-notation. */ /* -------------------------------------------------------------------- */ for (i = 0; i < n+1; i++) { ia[i] -= 1; } for (i = 0; i < nnz; i++) { ja[i] -= 1; } /* -------------------------------------------------------------------- */ /* .. 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, dparm); return 0; }