s_refine_gmres.c

Go to the documentation of this file.

/**
 *
 * @file s_refine_gmres.c
 *
 * PaStiX refinement functions implementations.
 *
 * @copyright 2015-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 *
 * @version 6.4.0
 * @author Mathieu Faverge
 * @author Pierre Ramet
 * @author Theophile Terraz
 * @author Xavier Lacoste
 * @author Gregoire Pichon
 * @author Vincent Bridonneau
 * @date 2024-07-05
 * @generated from /builds/2mk6rsew/0/solverstack/pastix/refinement/z_refine_gmres.c, normal z -> s, Tue Feb 25 14:36:04 2025
 *
 **/
#include "common.h"
#include "cblas.h"
#include "bcsc/bcsc.h"
#include "s_refine_functions.h"
 
/**
 *******************************************************************************
 *
 * @ingroup pastix_refine
 *
 * s_gmres_smp - Function computing GMRES iterative refinement.
 *
 *******************************************************************************
 *
 * @param[in] pastix_data
 *          The PaStiX data structure that describes the solver instance.
 *
 * @param[out] xp
 *          The solution vector.
 *
 * @param[in] bp
 *          The right hand side member (only one).
 *******************************************************************************
 *
 * @return Number of iterations
 *
 *******************************************************************************/
pastix_int_t
s_gmres_smp( pastix_data_t *pastix_data,
             pastix_rhs_t   xp,
             pastix_rhs_t   bp )
{
    struct s_solver     solver;
    Clock               refine_clk;
    float *x = (float*)(xp->b);
    float *b = (float*)(bp->b);
    float *gmHi, *gmH;
    float *gmVi, *gmV;
    float *gmWi, *gmW;
    float *gmcos, *gmsin;
    float *gmG;
    float *sgmWi   = NULL;
#if defined(PASTIX_DEBUG_GMRES)
    float *dbg_x, *dbg_r, *dbg_G;
#endif
    float  tmp;
    pastix_fixdbl_t     t0, t3;
    float              eps, resid, resid_b;
    float              norm, normb, normx;
    pastix_int_t        n, im, im1, itermax;
    pastix_int_t        i, j,  ldw, iters;
    int                 outflag, inflag;
    int                 savemem = 0;
    int                 precond = 1;
 
    memset( &solver, 0, sizeof(struct s_solver) );
    s_refine_init( &solver, pastix_data );
 
    /* if ( pastix_data->bcsc->mtxtype == PastixSymmetric ) { */
    /*     /\* Check if we need dot for non symmetric matrices (CEA patch) *\/ */
    /*     solver.dot = &bvec_sdot_seq; */
    /* } */
 
    /* Get the parameters */
    n       = pastix_data->bcsc->n;
    im      = pastix_data->iparm[IPARM_GMRES_IM];
    im1     = im + 1;
    itermax = pastix_data->iparm[IPARM_ITERMAX];
    eps     = pastix_data->dparm[DPARM_EPSILON_REFINEMENT];
    ldw     = n;
 
    if ( !(pastix_data->steps & STEP_NUMFACT) ) {
        precond = 0;
    }
 
    if ((!precond) || savemem ) {
        ldw = 0;
    }
 
    gmcos = (float *)solver.malloc(im  * sizeof(float));
    gmsin = (float *)solver.malloc(im  * sizeof(float));
    gmG   = (float *)solver.malloc(im1 * sizeof(float));
 
    /**
     * H stores the h_{i,j} elements of the upper hessenberg matrix H (See Alg. 9.5 p 270)
     * V stores the v_{i} vectors
     * W stores the M^{-1} v_{i} vectors to avoid the application of the
     *          preconditioner on the output result (See line 11 of Alg 9.5)
     *
     * If no preconditioner is applied, or the user wants to save memory, W
     * stores only temporarily one vector for the Ax product (ldw is set to 0 to
     * reuse the same vector at each iteration)
     */
    gmH = (float *)solver.malloc(im * im1 * sizeof(float));
    gmV = (float *)solver.malloc(n  * im1 * sizeof(float));
    if (precond && (!savemem) ) {
        gmW = (float *)solver.malloc(n * im  * sizeof(float));
    }
    else {
        gmW = (float *)solver.malloc(n       * sizeof(float));
    }
    memset( gmH, 0, im * im1 * sizeof(float) );
 
#if defined(PASTIX_DEBUG_GMRES)
    dbg_x = (float *)solver.malloc(n   * sizeof(float));
    dbg_r = (float *)solver.malloc(n   * sizeof(float));
    dbg_G = (float *)solver.malloc(im1 * sizeof(float));
    solver.copy( pastix_data, n, x, dbg_x );
#endif
 
    normb = solver.norm( pastix_data, n, b );
    if ( normb == 0. ) {
        normb = 1;
    }
    normx = solver.norm( pastix_data, n, x );
 
    /* Allocating a vector at half-precision, NULL pointer otherwise */
    if ( pastix_data->iparm[IPARM_MIXED] )
    {
        sgmWi = solver.malloc( n * sizeof(float) );
    }
 
    clockInit(refine_clk);
    clockStart(refine_clk);
 
    /**
     * Algorithm from Iterative Methods for Sparse Linear systems, Y. Saad, Second Ed. p267-273
     *
     * The version implemented is the Right preconditioned algorithm.
     */
    outflag = 1;
    iters = 0;
    while (outflag)
    {
        /* Initialize v_{0} and w_{0} */
        gmVi = gmV;
 
        /* Compute r0 = b - A * x */
        solver.copy( pastix_data, n, b, gmVi );
        if ( normx > 0. ) {
            solver.spmv( pastix_data, PastixNoTrans, -1., x, 1., gmVi );
        }
 
        /* Compute resid = ||r0||_f */
        resid = solver.norm( pastix_data, n, gmVi );
        resid_b = resid / normb;
 
        /* If residual is small enough, exit */
        if ( resid_b <= eps )
        {
            outflag = 0;
            break;
        }
 
        /* Compute v0 = r0 / resid */
        tmp = (float)( 1.0 / resid );
        solver.scal( pastix_data, n, tmp, gmVi );
 
        gmG[0] = (float)resid;
        inflag = 1;
        i = -1;
        gmHi = gmH - im1;
        gmWi = gmW - ldw;
 
        while( inflag )
        {
            clockStop( refine_clk );
            t0 = clockGet();
 
            i++;
 
            /* Set H and W pointers to the beginning of columns i */
            gmHi = gmHi + im1;
            gmWi = gmWi + ldw;
 
            /* Backup v_{i} into w_{i} for the end */
            solver.copy( pastix_data, n, gmVi, gmWi );
 
            /* Compute w_{i} = M^{-1} v_{i} */
            if ( precond ) {
                solver.spsv( pastix_data, gmWi, sgmWi );
            }
 
            /* v_{i+1} = A (M^{-1} v_{i}) = A w_{i} */
            gmVi += n;
            solver.spmv( pastix_data, PastixNoTrans, 1.0, gmWi, 0., gmVi );
 
            /* Classical Gram-Schmidt */
            for (j=0; j<=i; j++)
            {
                /* Compute h_{j,i} = < v_{i+1}, v_{j} > */
                gmHi[j] = solver.dot( pastix_data, n, gmVi, gmV + j * n );
 
                /* Compute v_{i+1} = v_{i+1} - h_{j,i} v_{j} */
                solver.axpy( pastix_data, n, -1. * gmHi[j],  gmV + j * n, gmVi );
            }
 
            /* Compute || v_{i+1} ||_f */
            norm = solver.norm( pastix_data, n, gmVi );
            gmHi[i+1] = norm;
 
            /* Compute v_{i+1} = v_{i+1} / h_{i+1,i} iff h_{i+1,i} is not too small */
            if ( norm > 1e-50 )
            {
                tmp = (float)(1.0 / norm);
                solver.scal( pastix_data, n, tmp, gmVi );
            }
 
            /* Apply the previous Givens rotation to the new column (should call LAPACKE_srot_work())*/
            for (j=0; j<i;j++)
            {
                /*
                 * h_{j,  i} = cos_j * h_{j,  i} +      sin_{j}  * h_{j+1, i}
                 * h_{j+1,i} = cos_j * h_{j+1,i} - (sin_{j}) * h_{j,   i}
                 */
                tmp = gmHi[j];
                gmHi[j]   = gmcos[j] * tmp       +      gmsin[j]  * gmHi[j+1];
                gmHi[j+1] = gmcos[j] * gmHi[j+1] - (gmsin[j]) * tmp;
            }
 
            /*
             * Compute the new Givens rotation (srotg)
             *
             * t = sqrtf( h_{i,i}^2 + h_{i+1,i}^2 )
             * cos = h_{i,i}   / t
             * sin = h_{i+1,i} / t
             */
            {
                tmp = sqrtf( gmHi[i]   * gmHi[i] +
                             gmHi[i+1] * gmHi[i+1] );
 
                if ( fabsf(tmp) <= eps ) {
                    tmp = (float)eps;
                }
                gmcos[i] = gmHi[i]   / tmp;
                gmsin[i] = gmHi[i+1] / tmp;
            }
 
            /* Update the residuals (See p. 168, eq 6.35) */
            gmG[i+1] = -gmsin[i] * gmG[i];
            gmG[i]   =  gmcos[i] * gmG[i];
 
            /* Apply the last Givens rotation */
            gmHi[i] = gmcos[i] * gmHi[i] + gmsin[i] * gmHi[i+1];
 
            /* (See p. 169, eq 6.42) */
            resid = fabsf( gmG[i+1] );
 
            resid_b = resid / normb;
            iters++;
            if ( (i+1 >= im) ||
                 (resid_b <= eps) ||
                 (iters >= itermax) )
            {
                inflag = 0;
            }
 
            clockStop((refine_clk));
            t3 = clockGet();
            if ( ( pastix_data->iparm[IPARM_VERBOSE] > PastixVerboseNot ) &&
                 ( pastix_data->procnum == 0 ) ) {
                solver.output_oneiter( t0, t3, resid_b, iters );
 
#if defined(PASTIX_DEBUG_GMRES)
                {
                    float normr2;
 
                    /* Compute y_m = H_m^{-1} g_m (See p. 169) */
                    memcpy( dbg_G, gmG, im1 * sizeof(float) );
                    cblas_strsv( CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit,
                                 i+1, gmH, im1, dbg_G, 1 );
 
                    solver.copy( pastix_data, n, b, dbg_r );
                    solver.copy( pastix_data, n, x, dbg_x );
 
                    /* Accumulate the current v_m */
                    solver.gemv( pastix_data, n, i+1, 1.0, (precond ? gmW : gmV), n, dbg_G, 1.0, dbg_x );
 
                    /* Compute b - Ax */
                    solver.spmv( pastix_data, PastixNoTrans, -1., dbg_x, 1., dbg_r );
 
                    normr2 = solver.norm( pastix_data, n, dbg_r );
                    fprintf(stdout, OUT_ITERREFINE_ERR, normr2 / normb );
                }
#endif
            }
        }
 
        /* Compute y_m = H_m^{-1} g_m (See p. 169) */
        cblas_strsv( CblasColMajor, CblasUpper, CblasNoTrans, CblasNonUnit,
                     i+1, gmH, im1, gmG, 1 );
 
        /**
         * Compute x_m = x_0 + M^{-1} V_m y_m
         *             = x_0 +        W_m y_m
         */
        if (precond && savemem) {
            /**
             * Since we saved memory, we do not have (M^{-1} V_m) stored,
             * thus we compute:
             *     w = V_m y_m
             *     w = M^{-1} (V_m y_m)
             *     x = x0 + (M^{-1} (V_m y_m))
             */
            solver.gemv( pastix_data, n, i+1, 1.0, gmV, n, gmG, 0., gmW );
            solver.spsv( pastix_data, gmW, sgmWi );
            solver.axpy( pastix_data, n, 1.,  gmW, x );
        }
        else {
            /**
             * Since we did not saved memory, we do have (M^{-1} V_m) stored in
             * W_m if precond is true, thus we compute:
             *     x = x0 + W_m y_m, if precond
             *     x = x0 + V_m y_m, if not precond
             */
            gmWi = precond ? gmW : gmV;
            solver.gemv( pastix_data, n, i+1, 1.0, gmWi, n, gmG, 1.0, x );
        }
 
        /**
         * Exit only if maximum number of iteration is reached.
         * Exit on residual if checked at the beginning of the outer loop to be
         * sure that the final residual of Ax-b is equal to the estimator
         * computed within the inner loop.
         */
        if (iters >= itermax)
        {
            outflag = 0;
        }
    }
 
    clockStop( refine_clk );
    t3 = clockGet();
 
    solver.output_final( pastix_data, resid_b, iters, t3, x, x );
 
    solver.free(gmcos);
    solver.free(gmsin);
    solver.free(gmG);
    solver.free(gmH);
    solver.free(gmV);
    solver.free(gmW);
    solver.free(sgmWi);
#if defined(PASTIX_DEBUG_GMRES)
    solver.free(dbg_x);
    solver.free(dbg_r);
    solver.free(dbg_G);
#endif
 
    return iters;
}