# Examples for Sparse Unsymmetric Linear Systems

In this section two examples (Fortran, C) are provided to solve unsymmetric linear systems with PARDISO. To solve the systems of equations `Ax = b`, where

## Example Results for Unsymmetric Systems

Upon successful execution of the solver, the result of the solution X is as follows

Reordering completed ...

Number of nonzeros in factors = 21

Number of factorization MFLOPS = 0

Factorization completed ...

Solve completed ...

The solution of the system is

x( 1) = -0.522321429

x( 2) = -0.00892857143

x( 3) = 1.22098214

x( 4) = -0.504464286

x( 5) = -0.214285714

## Example. pardiso_unsym.f for Unsymmetric Linear Systems

*******************************************************************************

* INTEL CONFIDENTIAL

* Copyright(C) 2004-2008 Intel Corporation. All Rights Reserved.

* The source code contained or described herein and all documents related to

* the source code ("Material") are owned by Intel Corporation or its suppliers

* or licensors. Title to the Material remains with Intel Corporation or its

* suppliers and licensors. The Material contains trade secrets and proprietary

* and confidential information of Intel or its suppliers and licensors. The

* Material is protected by worldwide copyright and trade secret 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, trade secret or other intellectual

* property right is granted to or conferred upon you by disclosure or delivery

* of the Materials, either expressly, by implication, inducement, estoppel or

* otherwise. Any license under such intellectual property rights must be

* express and approved by Intel in writing.

*

********************************************************************************

* Content : MKL PARDISO Fortran 77 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 http://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---------------------------------------------------------------------

PROGRAM pardiso_unsym

IMPLICIT NONE

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

INTEGER*8 pt(64)

C.. All other variables

INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl

INTEGER iparm(64)

INTEGER ia(6)

INTEGER ja(13)

REAL*8 a(13)

REAL*8 b(5)

REAL*8 x(5)

INTEGER i, idum

REAL*8 waltime1, waltime2, ddum

C.. Fill all arrays containing matrix data.

DATA n /5/, nrhs /1/, maxfct /1/, mnum /1/

DATA ia /1,4,6,9,12,14/

DATA ja

1 / 1, 2, 4,

2 1, 2,

3 3, 4, 5,

4 1, 3, 4,

5 2, 5/

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/

integer omp_get_max_threads

external omp_get_max_threads

C..

C.. Set up 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) = mkl_get_max_threads() ! numbers of processors, value of MKL_NUM_THREADS

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 compoments 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 ! perturbe the pivot elements with 1E-13

iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS

iparm(12) = 0 ! not in use

iparm(13) = 0 ! not in use

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.. Initiliaze the internal solver memory pointer. This is only

C necessary for the FIRST call of the PARDISO solver.

do i = 1, 64

pt(i) = 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,

1 idum, nrhs, iparm, msglvl, ddum, ddum, error)

WRITE(*,*) 'Reordering completed ... '

IF (error .NE. 0) THEN

WRITE(*,*) 'The following ERROR was detected: ', error

STOP

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,

1 idum, nrhs, iparm, msglvl, ddum, ddum, error)

WRITE(*,*) 'Factorization completed ... '

IF (error .NE. 0) THEN

WRITE(*,*) 'The following ERROR was detected: ', error

STOP

ENDIF

C.. Back substitution and iterative refinement

iparm(8) = 2 ! max numbers of iterative refinement steps

phase = 33 ! only factorization

do i = 1, n

b(i) = 1.d0

end do

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,

1 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,

1 idum, nrhs, iparm, msglvl, ddum, ddum, error)

END

## Example. pardiso_unsym.c for Unsymmetric Linear Systems

/*

********************************************************************************

* INTEL CONFIDENTIAL

* Copyright(C) 2004-2008 Intel Corporation. All Rights Reserved.

* The source code contained or described herein and all documents related to

* the source code ("Material") are owned by Intel Corporation or its suppliers

* or licensors. Title to the Material remains with Intel Corporation or its

* suppliers and licensors. The Material contains trade secrets and proprietary

* and confidential information of Intel or its suppliers and licensors. The

* Material is protected by worldwide copyright and trade secret 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, trade secret or other intellectual

* property right is granted to or conferred upon you by disclosure or delivery

* of the Materials, either expressly, by implication, inducement, estoppel or

* otherwise. Any license under such intellectual property rights must be

* express and approved by Intel in writing.

*

********************************************************************************

* 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: */

/* http://www.pardiso-project.org*/

/* */

/* (C) Olaf Schenk, Department of Computer Science, */

/* University of Basel, Switzerland.*/

/* Email: olaf.schenk@unibas.ch*/

/* --------------------------------------------------------------------*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

extern int omp_get_max_threads();

/* PARDISO prototype. */

#if defined(_WIN32) || defined(_WIN64)

#define pardiso_ PARDISO

#else

#define PARDISO pardiso_

#endif

extern int PARDISO

(void *, int *, int *, int *, int *, int *,

double *, int *, int*, int *, int *, int *,

int *, double *, double*, int *);

int main( void ) {

/* Matrix data. */

int n = 5;

int ia[ 6] = { 1, 4, 6, 9, 12, 14 };

int ja[13] = { 1, 2, 4,

1, 2,

3, 4, 5,

1, 3, 4,

2, 5 };

double a[18] = { 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 };

int mtype = 11; /* Real unsymmetric matrix */

/* RHS and solution vectors.*/

double b[5], x[5];

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];

int maxfct, mnum, phase, error, msglvl;

/* Auxiliary variables.*/

int i;

double ddum; /* Double dummy */

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 MKL_NUM_THREADS */

iparm[2] = mkl_get_max_threads();

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] = 0; /* Not in use */

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 = 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: %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");

/* --------------------------------------------------------------------*/

/* .. 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;

}

Parent topic: PARDISO Code Examples