/* ******************************************************************************** * Copyright(C) 2004-2011 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 symmetric 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" MKL_INT main( void ) { /* Matrix data. */ MKL_INT n = 27; MKL_INT ia[ 28] = { 1, 5, 10, 14, 19, 25, 30, 34, 39, 43, 48, 54, 59, 65, 72, 78, 83, 89, 94, 98, 103, 107, 112, 118, 123, 127, 132, 136 }; MKL_INT ja[137] = { 1, 2, 4, 10, 1, 2, 3, 5, 11, 2, 3, 6, 12, 1, 4, 5, 7, 13, 2, 4, 5, 6, 8, 14, 3, 5, 6, 9, 15, 4, 7, 8, 16, 5, 7, 8, 9, 17, 6, 8, 9, 18, 1, 10, 11, 13, 19, 2, 10, 11, 12, 14, 20, 3, 11, 12, 15, 21, 4, 10, 13, 14, 16, 22, 5, 11, 13, 14, 15, 17, 23, 6, 12, 14, 15, 18, 24, 7, 13, 16, 17, 25, 8, 14, 16, 17, 18, 26, 9, 15, 17, 18, 27, 10, 19, 20, 22, 11, 19, 20, 21, 23, 12, 20, 21, 24, 13, 19, 22, 23, 25, 14, 20, 22, 23, 24, 26, 15, 21, 23, 24, 27, 16, 22, 25, 26, 17, 23, 25, 26, 27, 18, 24, 26, 27 }; double a[137] = { 1.8199702E+17, -1.0985973E+16, -8.1912757E+14, -8.6409330E+15, -1.1592385E+16, 1.4588944E+17, -1.6598948E+14, -5.6100655E+14, -4.7898050E+15, -2.0017022E+14, 3.1136505E+15, -1.9205374E+13, -5.9795118E+13, -1.3057714E+15, 1.4310660E+16, -4.9571076E+14, -2.7139611E+11, -3.1083547E+14, -8.6506768E+14, -5.1366198E+14, 1.1123487E+16, -2.9914416E+13, -3.2975526E+11, -2.3459568E+14, -2.7025094E+13, -3.1866762E+13, 7.2196575E+14, -1.8579852E+11, -1.1541145E+13, -3.2787746E+11, 1.8697716E+13, -1.6858925E+12, -3.9060259E+11, -3.8703976E+11, -1.6760384E+12, 2.5464555E+13, -8.2631020E+11, -7.2688127E+11, -2.3000133E+11, -5.9658456E+11, 1.8608376E+13, -2.8582235E+11, -8.2937740E+15, 3.5702905E+17, -2.6536052E+16, -2.0436672E+15, -6.9282499E+15, -4.5801451E+15, -2.7880778E+16, 3.0043138E+17, -3.1712260E+14, -1.2732341E+15, -3.8033412E+15, -5.6776276E+13, -3.7229968E+14, 5.6106350E+15, -3.3883381E+13, -4.8613879E+13, -3.0381478E+14, -3.1016079E+15, 2.9225968E+16, -1.0906780E+15, -9.5587705E+10, -2.5677595E+14, -2.2928706E+14, -1.9219051E+15, -1.1272249E+15, 2.2035849E+16, -5.5743823E+13, -1.3566280E+11, -1.9474403E+14, -1.1203748E+13, -4.7868887E+13, -5.8888569E+13, 1.2579326E+15, -3.0559918E+11, -9.7358394E+12, -4.0293109E+11, -1.1471707E+11, 8.5987846E+12, -2.0592403E+12, -3.6577552E+11, -7.3020709E+11, -1.5613501E+11, -2.0062270E+12, 1.3477062E+13, -1.4252345E+12, -6.7511668E+11, -2.7899835E+11, -3.7727915E+11, -8.9305330E+11, 2.9937778E+13, -2.5662416E+11, -7.3386785E+15, 1.5517277E+17, -9.9426045E+15, -7.5088987E+14, -4.0517559E+15, -1.0494409E+16, 1.2732889E+17, -1.5332513E+14, -5.1321598E+14, -5.1979443E+13, -1.8264756E+14, 2.8940416E+15, -1.8055436E+13, -2.6498432E+14, -1.1933087E+15, 1.2956211E+16, -4.5460179E+14, -2.1820277E+11, -2.0089857E+14, -7.8911716E+14, -4.7076781E+14, 1.0125458E+16, -2.8321281E+13, -2.6898506E+11, -1.0108888E+13, -2.5336910E+13, -2.9912136E+13, 6.8683701E+14, -1.9152165E+11, -3.6779567E+11, -2.5767510E+11, 1.6941785E+13, -1.7538799E+12, -6.8959640E+11, -3.1031525E+11, -1.7352727E+12, 2.3639041E+13, -9.2134723E+11, -2.6378419E+11, -2.3387315E+11, -6.5677925E+11, 2.0127326E+13 }; MKL_INT mtype = 11; /* Real unsymmetric matrix */ /* RHS and solution vectors. */ double b[27], x[27]; 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. */ /* -------------------------------------------------------------------- */ /* .. 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 */ /* Numbers of processors, value of OMP_NUM_THREADS */ iparm[2] = 1; 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; /* Not in use */ 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 */ 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); 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); exit(2); } printf("\nFactorization completed ... ");*/ /* -------------------------------------------------------------------- */ /* .. Back substitution and iterative refinement. */ /* -------------------------------------------------------------------- */ phase = 23; iparm[3] = 0; // iparm[7] = 2; /* Max numbers of iterative refinement steps. */ /* Set right hand side to one. */ for (i = 0; i < n; i++) { b[i] = 1.E+15; } /* 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); exit(3); }*/ for (i = 0; i < 3; i++) { double a1[137] = { 3.1020985E+17, -2.4146940E+16, -2.1968032E+15, -1.0630049E+16, -1.4169095E+16, 2.2792357E+17, -3.9137316E+14, -1.4124705E+15, -4.2461123E+15, -2.5345015E+14, 4.6284779E+15, -5.1097062E+13, -4.0183219E+13, -9.6843386E+14, 2.1218966E+16, -9.7970278E+14, -7.7146096E+11, -2.5265583E+14, -6.0402454E+14, -5.5577625E+14, 1.5947305E+16, -6.6088637E+13, -8.8817710E+11, -1.8314715E+14, -1.6257566E+13, -3.4204870E+13, 1.0365083E+15, -7.0787052E+11, -6.9038679E+12, -1.8506283E+11, 3.0249444E+13, -3.6607025E+12, -8.1247935E+11, -2.1401515E+11, -2.8427986E+12, 4.0046829E+13, -1.8596013E+12, -1.3092627E+12, -1.8674792E+10, -9.4394509E+11, 2.7508076E+13, -1.8443366E+11, -2.3771729E+16, 6.3843459E+17, -6.0378852E+16, -5.3123017E+15, -5.9758409E+15, -1.2098926E+16, -3.7585413E+16, 4.8635701E+17, -7.7988922E+14, -3.1131311E+15, -2.2053970E+15, -1.5142300E+14, -5.2297777E+14, 8.5245190E+15, -9.2258297E+13, -1.7459978E+13, -7.7167167E+14, -2.5212207E+15, 4.3910479E+16, -2.1396915E+15, -2.0338260E+11, -1.2375781E+14, -5.5201258E+14, -1.3444418E+15, -1.2692116E+15, 3.1872486E+16, -1.2595582E+14, -2.9939794E+11, -9.1343578E+13, -2.9288806E+13, -2.7773192E+13, -6.7652374E+13, 1.8287424E+15, -1.1925013E+12, -2.7521088E+12, -1.2716220E+12, -3.0258489E+11, 1.6516554E+13, -5.0188832E+12, -5.5676056E+11, -1.9836806E+12, -3.5691227E+11, -4.6337475E+12, 2.4912287E+13, -3.2328951E+12, -9.1093657E+11, -7.5204375E+11, -3.6856676E+10, -1.8503116E+12, 4.5168795E+13, -7.9658402E+10, -1.9826735E+16, 2.8068455E+17, -2.3839661E+16, -2.1102090E+15, -1.0183662E+16, -1.4106058E+16, 2.0915494E+17, -3.8676744E+14, -1.3491130E+15, -1.3040333E+14, -2.5076093E+14, 4.4525413E+15, -5.0248332E+13, -6.2180423E+14, -7.8546805E+14, 1.9657415E+16, -9.3759753E+14, -6.0587455E+11, -4.4623777E+14, -4.8732175E+14, -5.3075370E+14, 1.4850881E+16, -6.5598134E+13, -7.1351199E+11, -2.4770874E+13, -1.3571612E+13, -3.2975480E+13, 1.0101377E+15, -7.5603522E+11, -1.0671994E+12, -2.1115614E+11, 2.8289268E+13, -3.9466851E+12, -1.7914308E+12, -2.3478239E+11, -3.1616399E+12, 3.8514771E+13, -2.1625491E+12, -6.7463289E+11, -3.0868685E+10, -1.1446616E+12, 3.0518366E+13 }; if (i>0) { iparm[3] = 31; /* No iterative-direct algorithm */ PARDISO (pt, &maxfct, &mnum, &mtype, &phase, &n, a1, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error); } else { 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); exit(3); } printf("\nSolve completed ... "); printf("\nThe solution of the system is: "); for (j = 0; j < n; j++) { printf("\n x [%d] = % f", j, x[j] ); } 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; }