PaStiX Handbook  6.2.1
schur.c File Reference

Schur usage example. More...

Go to the source code of this file.

Detailed Description

Schur usage example.

Version
6.2.0
Author
Pierre Ramet
Mathieu Faverge
Matias Hastaran
Tony Delarue
Date
2021-04-07
/
#include <pastix.h>
#include <spm.h>
#include <lapacke.h>
void
schurFactorize( pastix_coeftype_t flttype,
pastix_factotype_t factotype,
pastix_int_t N,
void *S,
pastix_int_t lds,
int **ipiv )
{
int info = 0;
assert( ipiv != NULL );
if ( factotype == PastixFactGETRF ) {
*ipiv = malloc( N * sizeof(int) );
}
switch (flttype) {
case PastixFloat:
switch (factotype) {
info = LAPACKE_spotrf_work( LAPACK_COL_MAJOR, 'L', N, S, lds );
break;
info = LAPACKE_sgetrf_work( LAPACK_COL_MAJOR, N, N, S, lds, *ipiv );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
break;
case PastixComplex32:
switch (factotype) {
info = LAPACKE_cpotrf_work( LAPACK_COL_MAJOR, 'L', N, S, lds );
break;
info = LAPACKE_cgetrf_work( LAPACK_COL_MAJOR, N, N, S, lds, *ipiv );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
break;
case PastixComplex64:
switch (factotype) {
info = LAPACKE_zpotrf_work( LAPACK_COL_MAJOR, 'L', N, S, lds );
break;
info = LAPACKE_zgetrf_work( LAPACK_COL_MAJOR, N, N, S, lds, *ipiv );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
break;
case PastixDouble:
switch (factotype) {
info = LAPACKE_dpotrf_work( LAPACK_COL_MAJOR, 'L', N, S, lds );
break;
info = LAPACKE_dgetrf_work( LAPACK_COL_MAJOR, N, N, S, lds, *ipiv );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
break;
default:
fprintf(stderr, "Incorrect arithmetic type\n");
}
if (info != 0) {
fprintf(stderr, "Error in schurFactorize with info =%d\n", info );
}
return;
}
void
schurSolve( pastix_coeftype_t flttype,
pastix_factotype_t factotype,
pastix_int_t N,
pastix_int_t Nschur,
pastix_int_t NRHS,
void *S,
pastix_int_t lds,
void *bptr,
pastix_int_t ldb,
int **ipiv )
{
int info = 0;
assert(ipiv != NULL);
switch (flttype) {
case PastixFloat:
{
float *b = (float *)bptr;
b += N - Nschur;
switch (factotype) {
info = LAPACKE_spotrs_work( LAPACK_COL_MAJOR, 'L', Nschur, NRHS, S, lds, b, ldb );
break;
info = LAPACKE_sgetrs_work( LAPACK_COL_MAJOR, 'N', Nschur, NRHS, S, lds, *ipiv, b, ldb );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
}
break;
case PastixComplex32:
{
pastix_complex32_t *b = (pastix_complex32_t *)bptr;
b += N - Nschur;
switch (factotype) {
info = LAPACKE_cpotrs_work( LAPACK_COL_MAJOR, 'L', Nschur, NRHS, S, lds, b, ldb );
break;
info = LAPACKE_cgetrs_work( LAPACK_COL_MAJOR, 'N', Nschur, NRHS, S, lds, *ipiv, b, ldb );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
}
break;
case PastixComplex64:
{
pastix_complex64_t *b = (pastix_complex64_t *)bptr;
b += N - Nschur;
switch (factotype) {
info = LAPACKE_zpotrs_work( LAPACK_COL_MAJOR, 'L', Nschur, NRHS, S, lds, b, ldb );
break;
info = LAPACKE_zgetrs_work( LAPACK_COL_MAJOR, 'N', Nschur, NRHS, S, lds, *ipiv, b, ldb );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
}
break;
case PastixDouble:
{
double *b = (double *)bptr;
b += N - Nschur;
switch (factotype) {
info = LAPACKE_dpotrs_work( LAPACK_COL_MAJOR, 'L', Nschur, NRHS, S, lds, b, ldb );
break;
info = LAPACKE_dgetrs_work( LAPACK_COL_MAJOR, 'N', Nschur, NRHS, S, lds, *ipiv, b, ldb );
break;
default:
fprintf(stderr, "Factorization type not handled by Schur example\n");
}
}
break;
default:
fprintf(stderr, "Incorrect arithmetic type\n");
}
if (*ipiv != NULL) {
free( *ipiv );
*ipiv = NULL;
}
if (info != 0) {
fprintf(stderr, "Error in schurSolve with info =%d\n", info );
}
return;
}
int main (int argc, char **argv)
{
pastix_data_t *pastix_data = NULL; /*< Pointer to the storage structure required by pastix */
pastix_int_t iparm[IPARM_SIZE]; /*< Integer in/out parameters for pastix */
double dparm[DPARM_SIZE]; /*< Floating in/out parameters for pastix */
spm_driver_t driver;
char *filename;
spmatrix_t *spm, spm2;
void *x, *b, *S, *x0 = NULL;
size_t size;
int check = 1;
int nrhs = 1;
int rc = 0;
pastix_int_t nschur, lds, ldb;
int *ipiv = NULL;
/**
* Initialize parameters to default values
*/
pastixInitParam( iparm, dparm );
/**
* Get options from command line
*/
pastixGetOptions( argc, argv,
iparm, dparm,
&check, &driver, &filename );
if ( (iparm[IPARM_FACTORIZATION] == PastixFactLDLT) ||
{
fprintf(stderr, "This types of factorization (LDL^t and LDL^h) are not supported by this example.\n");
return EXIT_FAILURE;
}
/**
* Startup PaStiX
*/
pastixInit( &pastix_data, MPI_COMM_WORLD, iparm, dparm );
/**
* Read the sparse matrix with the driver
*/
spm = malloc( sizeof( spmatrix_t ) );
spmReadDriver( driver, filename, spm );
free( filename );
spmPrintInfo( spm, stdout );
rc = spmCheckAndCorrect( spm, &spm2 );
if ( rc != 0 ) {
spmExit( spm );
*spm = spm2;
rc = 0;
}
/**
* Generate a Fake values array if needed for the numerical part
*/
if ( spm->flttype == SpmPattern ) {
spmGenFakeValues( spm );
}
/**
* Initialize the schur list with the first third of the unknowns
*/
{
nschur = spm->gN / 3;
/* Set to a maximum to avoid memory problem with the test */
nschur = (nschur > 5000) ? 5000 : nschur;
if ( nschur > 0 ) {
pastix_int_t i;
pastix_int_t baseval = spmFindBase(spm);
pastix_int_t *list = (pastix_int_t*)malloc(nschur * sizeof(pastix_int_t));
for (i=0; i<nschur; i++) {
list[i] = i+baseval;
}
pastixSetSchurUnknownList( pastix_data, nschur, list );
free( list );
}
iparm[IPARM_SCHUR_SOLV_MODE] = PastixSolvModeInterface;
}
/**
* Perform ordering, symbolic factorization, and analyze steps
*/
pastix_task_analyze( pastix_data, spm );
/**
* Normalize A matrix (optional, but recommended for low-rank functionality)
*/
double normA = spmNorm( SpmFrobeniusNorm, spm );
spmScalMatrix( 1./normA, spm );
/**
* Perform the numerical factorization
*/
pastix_task_numfact( pastix_data, spm );
/**
* Get the Schur complement back
*/
lds = nschur;
S = malloc( pastix_size_of( spm->flttype ) * nschur * lds );
rc = pastixGetSchur( pastix_data, S, lds );
if( rc == -1 ){
spmExit( spm );
free( spm );
free( S );
pastixFinalize( &pastix_data );
exit(0);
}
/**
* Factorize the Schur complement
*/
schurFactorize( spm->flttype, iparm[IPARM_FACTORIZATION],
nschur, S, lds, &ipiv );
/**
* Generates the b and x vector such that A * x = b
* Compute the norms of the initial vectors if checking purpose.
*/
size = pastix_size_of( spm->flttype ) * spm->n * nrhs;
x = malloc( size );
b = malloc( size );
ldb = spm->n;
if ( check )
{
if ( check > 1 ) {
x0 = malloc( size );
}
spmGenRHS( SpmRhsRndX, nrhs, spm, x0, spm->n, b, spm->n );
memcpy( x, b, size );
}
else {
spmGenRHS( SpmRhsRndB, nrhs, spm, NULL, spm->n, x, spm->n );
/* Apply also normalization to b vectors */
spmScalVector( spm->flttype, 1./normA, spm->n * nrhs, b, 1 );
}
/**
* Solve the linear system Ax = (P^tLUP)x = b
*/
/* 1- Apply P to b */
pastix_subtask_applyorder( pastix_data, spm->flttype,
PastixDirForward, spm->n, nrhs, x, ldb );
/* 2- Forward solve on the non Schur complement part of the system */
diag = PastixNonUnit;
}
else if( iparm[IPARM_FACTORIZATION] == PastixFactGETRF ) {
diag = PastixUnit;
}
pastix_subtask_trsm( pastix_data, spm->flttype,
nrhs, x, ldb );
/* 3- Solve the Schur complement part */
schurSolve( spm->flttype, iparm[IPARM_FACTORIZATION],
spm->n, nschur, nrhs, S, lds, x, ldb, &ipiv );
/* 4- Backward solve on the non Schur complement part of the system */
pastix_subtask_trsm( pastix_data, spm->flttype,
nrhs, x, ldb );
}
else if( iparm[IPARM_FACTORIZATION] == PastixFactGETRF ) {
pastix_subtask_trsm( pastix_data, spm->flttype,
nrhs, x, ldb );
}
/* 5- Apply P^t to x */
pastix_subtask_applyorder( pastix_data, spm->flttype,
PastixDirBackward, spm->n, nrhs, x, ldb );
if ( check )
{
rc = spmCheckAxb( dparm[DPARM_EPSILON_REFINEMENT], nrhs, spm, x0, spm->n, b, spm->n, x, spm->n );
if ( x0 ) {
free( x0 );
}
}
spmExit( spm );
free( spm );
free( S );
free( x );
free( b );
pastixFinalize( &pastix_data );
return rc;
}
/**
*

Definition in file schur.c.

pastix_subtask_applyorder
int pastix_subtask_applyorder(pastix_data_t *pastix_data, pastix_coeftype_t flttype, pastix_dir_t dir, pastix_int_t m, pastix_int_t n, void *b, pastix_int_t ldb)
Apply a permutation on the right-and-side vector before the solve step.
Definition: pastix_task_solve.c:122
IPARM_SCHUR_SOLV_MODE
@ IPARM_SCHUR_SOLV_MODE
Definition: api.h:106
DPARM_EPSILON_REFINEMENT
@ DPARM_EPSILON_REFINEMENT
Definition: api.h:154
pastixSetSchurUnknownList
void pastixSetSchurUnknownList(pastix_data_t *pastix_data, pastix_int_t n, const pastix_int_t *list)
Set the list of unknowns that belongs to the schur complement.
Definition: schur.c:49
pastixGetOptions
void pastixGetOptions(int argc, char **argv, pastix_int_t *iparm, double *dparm, int *check, spm_driver_t *driver, char **filename)
PaStiX helper function to read command line options in examples.
Definition: get_options.c:142
PastixLower
@ PastixLower
Definition: api.h:443
pastix.h
PastixDirForward
@ PastixDirForward
Definition: api.h:489
PastixNoTrans
@ PastixNoTrans
Definition: api.h:422
PastixUpper
@ PastixUpper
Definition: api.h:442
PastixConjTrans
@ PastixConjTrans
Definition: api.h:424
PastixNonUnit
@ PastixNonUnit
Definition: api.h:463
PastixFactGETRF
@ PastixFactGETRF
Definition: api.h:297
pastix_factotype_t
enum pastix_factotype_e pastix_factotype_t
Factorization algorithms available for IPARM_FACTORIZATION parameter.
PastixDirBackward
@ PastixDirBackward
Definition: api.h:490
pastixInit
void pastixInit(pastix_data_t **pastix_data, PASTIX_Comm pastix_comm, pastix_int_t *iparm, double *dparm)
Initialize the solver instance.
Definition: api.c:857
PastixUnit
@ PastixUnit
Definition: api.h:464
PastixFactLDLT
@ PastixFactLDLT
Definition: api.h:301
PastixFactPOTRF
@ PastixFactPOTRF
Definition: api.h:295
pastix_diag_t
enum pastix_diag_e pastix_diag_t
Diagonal.
main
int main(int argc, char *argv[])
Definition: binding_for_multimpi.c:49
pastixInitParam
void pastixInitParam(pastix_int_t *iparm, double *dparm)
Initialize the iparm and dparm arrays to their default values.
Definition: api.c:401
pastixGetSchur
int pastixGetSchur(const pastix_data_t *pastix_data, void *S, pastix_int_t lds)
Return the Schur complement.
Definition: schur.c:89
pastix_task_analyze
int pastix_task_analyze(pastix_data_t *pastix_data, const spmatrix_t *spm)
Perform all the preprocessing steps: ordering, symbolic factorization, reordering,...
Definition: pastix_task_analyze.c:53
pastixFinalize
void pastixFinalize(pastix_data_t **pastix_data)
Finalize the solver instance.
Definition: api.c:880
pastix_coeftype_t
#define pastix_coeftype_t
Arithmetic types.
Definition: api.h:281
pastix_subtask_trsm
int pastix_subtask_trsm(pastix_data_t *pastix_data, pastix_coeftype_t flttype, pastix_side_t side, pastix_uplo_t uplo, pastix_trans_t trans, pastix_diag_t diag, pastix_int_t nrhs, void *b, pastix_int_t ldb)
Apply a triangular solve on the right-and-side vectors.
Definition: pastix_task_solve.c:236
PastixLeft
@ PastixLeft
Definition: api.h:471
pastix_task_numfact
int pastix_task_numfact(pastix_data_t *pastix_data, spmatrix_t *spm)
Perform all the numerical factorization steps: fill the internal block CSC and the solver matrix stru...
Definition: pastix_task_sopalin.c:587
PastixFactLDLH
@ PastixFactLDLH
Definition: api.h:304
IPARM_FACTORIZATION
@ IPARM_FACTORIZATION
Definition: api.h:99