Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
6434 Discussions

Schur complement for asymmetric matrix

Hainan_W
Beginner
121 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
121 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

 

Gennady_F_Intel
Moderator
121 Views

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

Hainan_W
Beginner
121 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

Gennady_F_Intel
Moderator
121 Views

the same results. 

Reply