Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Schur complement for asymmetric matrix

Hainan_W
Beginner
593 Views

Hi all,

I am using mkl pardiso to do Schur complement. So far I got it working for symmetric matrix.  But for asymmetric matrix, I got crash at phase = 11. The error message is "... Access violation writing location ...." 

Can anyone provide any guidance on what the possible root cause is? I am using mkl 2018 update 2. The source code is as below. 

Thanks!

Hainan

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mkl.h>
#include <iostream>
 
void main()
{
int n = 5;
double a[11] = { 2, 1.0,
1.0, 1.0,
6.0, 7.0,
2.5,
1.0, 1.0, 7.0, 3.5 };
int ia[6] = { 0, 2, 4, 6, 7, 11 };
int ja[11] = { 0, 4,
1, 4,
2, 4,
3,
0, 1, 2, 4 };
int matrix_order = LAPACK_ROW_MAJOR;
char uplo = 'U';
int mtype = 11;   /* Real asymmetric matrix */
  /* RHS and solution vectors. */
//double b[8], x[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];
int maxfct, mnum, phase, error, msglvl, info;
/* Auxiliary variables. */
int i, j;
double ddum;          /* Double dummy */
int idum;         /* Integer dummy. */
  /* Schur data */
double schur[3] = { 0.0, 0.0,
0.0 };
int perm[5] = { 0, 0, 1, 1, 1 };
int ipiv[2];
int n_schur = 3; /* Schur complement solution size */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
for (i = 0; i < 64; i++)
{
iparm = 0;
}
iparm[1 - 1] = 1;         /* No solver default */
iparm[2 - 1] = 2;         /* Fill-in reordering from METIS */
iparm[10 - 1] = 8;        /* Perturb the pivot elements with 1E-13 */
iparm[11 - 1] = 0;        /* Use nonsymmetric permutation and scaling MPS */
iparm[13 - 1] = 0;        /* Maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm[12] = 1 in case of inappropriate accuracy */
iparm[14 - 1] = 0;        /* Output: Number of perturbed pivots */
iparm[18 - 1] = -1;       /* Output: Number of nonzeros in the factor LU */
iparm[19 - 1] = -1;       /* Output: Mflops for LU factorization */
iparm[35 - 1] = 1;        /*0: one-base; 1: zero-base*/
iparm[36 - 1] = 1;        /* Use Schur complement */
 
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 = 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, perm, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
if (error != 0)
{
printf("\nERROR during symbolic factorization: %d", error);
exit(1);
}
}
0 Kudos
4 Replies
Gennady_F_Intel
Moderator
593 Views

..\MKL_Forums\u799124>test.exe

=== PARDISO: solving a real nonsymmetric system ===
0-based array is turned ON
PARDISO double precision computation is turned ON
User provided fill-in reducing permutation is turned ON


Summary: ( reordering phase )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000021 s
Time spent in reordering of the initial matrix (reorder)         : 0.000007 s
Time spent in symbolic factorization (symbfct)                   : 0.014133 s
Time spent in data preparations for factorization (parlist)      : 0.000002 s
Time spent in allocation of internal data structures (malloc)    : 0.022789 s
Time spent in additional calculations                            : 0.000098 s
Total time spent                                                 : 0.037050 s

Statistics:
===========
Parallel Direct Factorization is running on 2 OpenMP

< Linear system Ax = b >
             number of equations:           5
             number of non-zeros in A:      11
             number of non-zeros in A (%): 44.000000

             number of right-hand sides:    1

< Factors L and U >
             number of columns for each panel: 128
             number of independent subgraphs:  0
< Preprocessing with input permutation >
             number of supernodes:                    5
             size of largest supernode:               1
             number of non-zeros in L:                8
             number of non-zeros in U:                3
             number of non-zeros in L+U:              11

 

0 Kudos
Gennady_F_Intel
Moderator
593 Views

I compiled ( icl /Qmkl test.cpp) you case as is. it works. see the output see above. mkl v.2018 u4, win10. 

0 Kudos
Hainan_W
Beginner
593 Views

Hi Gennady,

Thanks for the test and update. I am using mkl 2018 u2, win7. Looks it crashes only in debug mode. Can you have a try in debug mode?

Regards,

Hainan

0 Kudos
Gennady_F_Intel
Moderator
593 Views

the same results. 

0 Kudos
Reply