/*******************************************************************************
* Copyright 2004-2016 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.
*
* 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 usage of the "PARDISO" routine */
/* for Schur complement solve on symmetric linear systems */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mkl.h"

MKL_INT main (void)
{
    /* Matrix data. */
    MKL_INT n = 6;
    MKL_INT ia[7] = {0, 2, 4, 5, 8, 9, 10};
    MKL_INT ja[10] = {0, 1, 1, 2, 2, 3, 4, 5, 4, 5};
    double a[10] = { 1.0, 2.0, 2.0, 4.0, 1.0, 2.0, 3.0, 5.0, 1.0, 1.0};
    char uplo = 'U';
    MKL_INT mtype = -2;   /* Real symmetric matrix */
    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, info;
    /* Auxiliary variables. */
    MKL_INT i, j;
    double ddum = 0;          /* Double dummy */
    /* Schur data */
    double schur[4] = {0.0, 0.0,
                       0.0, 0.0 };
    MKL_INT perm[6] = {0, 1, 0, 1, 0, 0};
    MKL_INT n_schur = 2; /* Schur complement solution size */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
  
	iparm[0] = 0;
	pardisoinit(pt, &mtype, iparm);
	
	iparm[34] = 1;	/* PARDISO use C-style indexing for columns and pRow arrays */
	iparm[35] = 1; /* Use Schur complement */
	
	nrhs = 1;
    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 */

/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates   */
/* all memory that is necessary for the factorization.                  */
/* -------------------------------------------------------------------- */
    int idum = 0;
	phase = 11;
    pardiso (pt, &maxfct, &mnum, &mtype, &phase,
             &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);	 
			 
    if ( error != 0 )
    {
        printf ("\nERROR during symbolic factorization: %d", error);
        exit (1);
    }
    printf ("\nReordering completed ... ");
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
	idum = 0;
    phase = 22;
    pardiso(pt, &maxfct, &mnum, &mtype, &phase,
            &n, a, ia, ja, perm, &nrhs,
            iparm, &msglvl, &ddum, &schur, &error);
			
    if ( error != 0 )
    {
        printf ("\nERROR during numerical factorization: %d", error);
        exit (2);
    }
    printf ("\nFactorization completed ... ");
    printf("\nSchur matrix: \n"); fflush(0);
    for(i=0; i<n_schur; i++)
    {
        for(j=0; j<n_schur; j++)
        {
            printf("%e ",schur[i*n_schur+j]); fflush(0);
        }
        printf("\n"); fflush(0);
    }
/* -------------------------------------------------------------------- */
/* .. 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;
}
