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 = NULL;
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 ) );
rc = spmReadDriver( driver, filename, spm );
free( filename );
if ( rc != SPM_SUCCESS ) {
pastixFinalize( &pastix_data );
return rc;
}
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 );
spmScal( 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->nexp * nrhs;
x = malloc( size );
b = malloc( size );
ldb = spm->nexp;
if ( check )
{
if ( check > 1 ) {
x0 = malloc( size );
}
spmGenRHS( SpmRhsRndX, nrhs, spm, x0, spm->nexp, b, spm->nexp );
memcpy( x, b, size );
}
else {
spmGenRHS( SpmRhsRndB, nrhs, spm, NULL, spm->nexp, x, spm->nexp );
/* Apply also normalization to b vectors */
spmScalMat( 1./normA, spm, nrhs, b, spm->nexp );
}
/**
* Solve the linear system Ax = (P^tLUP)x = b
*/
/* 1- Apply P to b */
pastix_subtask_applyorder( pastix_data, spm->flttype,
PastixDirForward, spm->nexp, 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->nexp, 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->nexp, nrhs, x, ldb );
if ( check )
{
rc = spmCheckAxb( dparm[DPARM_EPSILON_REFINEMENT], nrhs, spm, x0, spm->nexp, b, spm->nexp, x, spm->nexp );
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:445
pastix.h
PastixDirForward
@ PastixDirForward
Definition: api.h:491
PastixNoTrans
@ PastixNoTrans
Definition: api.h:424
PastixUpper
@ PastixUpper
Definition: api.h:444
PastixConjTrans
@ PastixConjTrans
Definition: api.h:426
PastixNonUnit
@ PastixNonUnit
Definition: api.h:465
PastixFactGETRF
@ PastixFactGETRF
Definition: api.h:299
pastix_factotype_t
enum pastix_factotype_e pastix_factotype_t
Factorization algorithms available for IPARM_FACTORIZATION parameter.
PastixDirBackward
@ PastixDirBackward
Definition: api.h:492
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:859
PastixUnit
@ PastixUnit
Definition: api.h:466
PastixFactLDLT
@ PastixFactLDLT
Definition: api.h:303
PastixFactPOTRF
@ PastixFactPOTRF
Definition: api.h:297
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:882
pastix_coeftype_t
#define pastix_coeftype_t
Arithmetic types.
Definition: api.h:283
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:473
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:623
PastixFactLDLH
@ PastixFactLDLH
Definition: api.h:306
IPARM_FACTORIZATION
@ IPARM_FACTORIZATION
Definition: api.h:99