/* ******************************************************************************** * Copyright(C) 2004-2015 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. * * *Third Party trademarks are the property of their respective owners. * * 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 C example * ******************************************************************************** */ /* -------------------------------------------------------------------- */ /* Example program to show the use of the "PARDISO" routine */ /* on unsymmetric linear systems */ /* -------------------------------------------------------------------- */ /* This program can be downloaded from the following site: */ /* www.pardiso-project.org */ /* */ /* (C) Olaf Schenk, Department of Computer Science, */ /* University of Basel, Switzerland. */ /* Email: olaf.schenk@unibas.ch */ /* -------------------------------------------------------------------- */ #include #include #include #include "mkl_pardiso.h" #include "mkl_types.h" #include "mkl_spblas.h" #define USE_ONE_INDEX 1 MKL_INT main (void) { /* Matrix data. */ MKL_INT n = 5; #if USE_ONE_INDEX MKL_INT ia[6] = { 1, 4, 6, 9, 12, 14}; MKL_INT ja[13] = { 1, 2, 4, 1, 2, 3, 4, 5, 1, 3, 4, 2, 5 }; #else MKL_INT ia[6] = { 0, 3, 5, 8, 11, 13}; MKL_INT ja[13] = { 0, 1, 3, 0, 1, 2, 3, 4, 0, 2, 3, 1, 4 }; #endif double a[13] = { 1.0,-1.0, -3.0, -2.0, 5.0, 4.0, 6.0, 4.0, -4.0, 2.0, 7.0, 8.0, -5.0 }; MKL_INT mtype = 11; /* Real unsymmetric matrix */ /* RHS and solution vectors. */ double b[5], x[5], bs[5], res, res0; 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, j; double ddum; /* Double dummy */ MKL_INT idum; /* Integer dummy. */ char *uplo; /* -------------------------------------------------------------------- */ /* .. 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[6] = 0; /* Not in use */ iparm[7] = 2; /* Max numbers of iterative refinement steps */ iparm[8] = 0; /* Not in use */ iparm[9] = 13; /* Perturb the pivot elements with 1E-13 */ iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ iparm[11] = 0; /* Conjugate transposed/transpose solve */ iparm[12] = 1; /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */ iparm[13] = 0; /* Output: Number of perturbed pivots */ iparm[14] = 0; /* Not in use */ iparm[15] = 0; /* Not in use */ iparm[16] = 0; /* Not in use */ 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 */ #if USE_ONE_INDEX #else iparm[34] = 1; #endif 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: %d", error); getchar(); exit (1); } printf ("\nReordering completed ... "); printf ("\nNumber of nonzeros in factors = %d", iparm[17]); printf ("\nNumber of factorization MFLOPS = %d", 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: %d", error); getchar(); exit (2); } printf ("\nFactorization completed ... "); /* -------------------------------------------------------------------- */ /* .. Back substitution and iterative refinement. */ /* -------------------------------------------------------------------- */ phase = 33; /* Set right hand side to one. */ for ( i = 0; i < n; i++ ) { b[i] = 1; } // Loop over 3 solving steps: Ax=b, AHx=b and ATx=b for ( i = 0; i < 3; i++ ) { iparm[11] = i; /* Conjugate transposed/transpose solve */ if ( i == 0 ) uplo = "non-transposed"; else if ( i == 1 ) uplo = "conjugate transposed"; else uplo = "transposed"; printf ("\n\nSolving %s system...\n", uplo); PARDISO (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error); if ( error != 0 ) { printf ("\nERROR during solution: %d", error); getchar(); exit (3); } printf ("\nThe solution of the system is: "); for ( j = 0; j < n; j++ ) { printf ("\n x [%d] = % f", j, x[j]); } printf ("\n"); // Compute residual mkl_dcsrgemv (uplo, &n, a, ia, ja, x, bs); res = 0.0; res0 = 0.0; for ( j = 1; j <= n; j++ ) { res += (bs[j - 1] - b[j - 1]) * (bs[j - 1] - b[j - 1]); res0 += b[j - 1] * b[j - 1]; } res = sqrt (res) / sqrt (res0); printf ("\nRelative residual = %e", res); // Check residual if ( res > 1e-10 ) { printf ("Error: residual is too high!\n"); getchar(); exit (10 + i); } } /* -------------------------------------------------------------------- */ /* .. 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); getchar(); return 0; }