Schur usage example.
/
#include <spm.h>
#include <lapacke.h>
void
pastix_int_t N,
void *S,
pastix_int_t lds,
int **ipiv )
{
int info = 0;
assert( ipiv != NULL );
*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
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;
pastix_int_t iparm[IPARM_SIZE];
double dparm[DPARM_SIZE];
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;
iparm, dparm,
&check, &driver, &filename );
{
fprintf(stderr, "This types of factorization (LDL^t and LDL^h) are not supported by this example.\n");
return EXIT_FAILURE;
}
pastixInit( &pastix_data, MPI_COMM_WORLD, iparm, dparm );
spm = malloc( sizeof( spmatrix_t ) );
rc = spmReadDriver( driver, filename, spm );
free( filename );
if ( rc != SPM_SUCCESS ) {
return rc;
}
spmPrintInfo( spm, stdout );
rc = spmCheckAndCorrect( spm, &spm2 );
if ( rc != 0 ) {
spmExit( spm );
*spm = spm2;
rc = 0;
}
if ( spm->flttype == SpmPattern ) {
spmGenFakeValues( spm );
}
{
nschur = spm->gN / 3;
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;
}
free( list );
}
}
double normA = spmNorm( SpmFrobeniusNorm, spm );
spmScal( 1./normA, spm );
lds = nschur;
S = malloc( pastix_size_of( spm->flttype ) * nschur * lds );
if( rc == -1 ){
spmExit( spm );
free( spm );
free( S );
exit(0);
}
nschur, S, lds, &ipiv );
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 );
spmScalMat( 1./normA, spm, nrhs, b, spm->nexp );
}
}
}
nrhs, x, ldb );
spm->nexp, nschur, nrhs, S, lds, x, ldb, &ipiv );
nrhs, x, ldb );
}
nrhs, x, ldb );
}
if ( check )
{
if ( x0 ) {
free( x0 );
}
}
spmExit( spm );
free( spm );
free( S );
free( x );
free( b );
return rc;
}