step-by-step.c File Reference

A step-by-step example that runs one full analyze (ordering, symbolic factorization, analyze), then loops over 2 factorizations that are both used for 2 solves each. More...

Go to the source code of this file.

Detailed Description

A step-by-step example that runs one full analyze (ordering, symbolic factorization, analyze), then loops over 2 factorizations that are both used for 2 solves each.

Copyright

2015-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria, Univ. Bordeaux. All rights reserved.

Version

6.4.0

Author

Mathieu Faverge

Pierre Ramet

Tony Delarue

Alycia Lisito

Date

2024-07-05

  /
#include <pastix.h>
#include <spm.h>
 
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, *x0 = NULL;
    size_t          size;
    int             scatter = 0;
    int             check   = 1;
    int             nrhs    = 10;
    int             rc      = PASTIX_SUCCESS;
    int             nfact   = 2;
    int             nsolv   = 2;
    long            i, j;
 
    /**
       Initialize parameters to default values
      /
    pastixInitParam( iparm, dparm );
 
    /**
       Get options from command line
      /
    pastixGetOptions( argc, argv,
                      iparm, dparm,
                      &check, &scatter, &driver, &filename );
 
    /**
       Startup PaStiX
      /
    pastixInit( &pastix_data, MPI_COMM_WORLD, iparm, dparm );
 
    /**
       Read the sparse matrix with the driver
      /
    spm = malloc( sizeof( spmatrix_t ) );
    if ( scatter ) {
        rc = spmReadDriverDist( driver, filename, spm, MPI_COMM_WORLD );
    }
    else {
        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 );
    }
 
    /**
       Perform ordering, symbolic factorization, and analyze steps
      /
    pastix_subtask_order( pastix_data, spm, NULL );
    pastix_subtask_symbfact( pastix_data );
    pastix_subtask_reordering( pastix_data );
    pastix_subtask_blend( pastix_data );
 
    /**
       Normalize A matrix (optional, but recommended for low-rank functionality)
      /
    double normA = spmNorm( SpmFrobeniusNorm, spm );
    spmScal( 1./normA, spm );
 
    size = pastix_size_of( spm->flttype ) * spm->nexp * nrhs;
    x = malloc( size );
    b = malloc( size );
    if ( check > 1 ) {
        x0 = malloc( size );
    }
 
    /* Do nfact factorization */
    for (i = 0; i < nfact; i++)
    {
        /**
           Perform the numerical factorization
          /
        pastix_subtask_spm2bcsc( pastix_data, spm );
        pastix_subtask_bcsc2ctab( pastix_data );
        pastix_subtask_sopalin( pastix_data );
        for (j = 0; j < nsolv; j++)
        {
            /**
               Generates the b and x vector such that A * x = b
               Compute the norms of the initial vectors if checking purpose.
              /
            if ( check )
            {
                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 );
 
                /* Save b for refinement */
                memcpy( b, x, size );
            }
 
            /**
               Solve the linear system
              /
            pastix_task_solve( pastix_data, spm->nexp, nrhs, x, spm->nexp );
            pastix_task_refine( pastix_data, spm->nexp, nrhs, b, spm->nexp, x, spm->nexp );
 
            if ( check ) {
                rc |= spmCheckAxb( dparm[DPARM_EPSILON_REFINEMENT], nrhs, spm, x0, spm->nexp, b, spm->nexp, x, spm->nexp );
            }
        }
    }
 
    spmExit( spm );
    free( spm );
    free( b );
    free( x );
    if ( x0 ) {
        free( x0 );
    }
    pastixFinalize( &pastix_data );
 
    return rc;
}
 
/**

Definition in file step-by-step.c.