krylov.c

Go to the documentation of this file.

 
#include "krylov.h"
#include "precon.h"
#include <math.h>
#include <omp.h>
 
 
ucfd_status_t serial_gmres(sparse_matrix_t op, ucfd_precon_type_t precon_type, int bn, int block, int m, int *iter, double tol,
                           int *row_ptr, int *col_ind, int *diag_ind, double *precon_nnz_data,
                           double *x, double *b, double *H, double *V, double *g, double *y, double *w, double *r)
{
    int it, i, j, itmax;
    const int n = bn * BLOCK;
    double beta, tmp, c, s, h1, h2, rr;
    ucfd_precon_solve psolve;
    sparse_status_t mklstat;
 
    /* --------------------------------
     * Set variables
    -------------------------------- */
    // Sparse matrix description for system matrix A(op)
    struct matrix_descr descr;
    descr.type = SPARSE_MATRIX_TYPE_GENERAL;
    descr.mode = 0;
    descr.diag = 0;
    
    // Set preconditioner solver
    if (precon_type == BILU) psolve = bilu_psolve;
    else if (precon_type == LUSGS) psolve = lusgs_psolve;
    else psolve = none_psolve;
 
    // Maximum iteration
    itmax = *iter;
    it = 0;
 
    /* --------------------------------
     * Initial residual
     1) r := -A @ x
     2) r := b -A @ x (r += rhs)
    -------------------------------- */
    mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1.0, op, descr, x, 0.0, r);
    if (mklstat != SPARSE_STATUS_SUCCESS) {
        *iter = mklstat;
        return UCFD_MKL_FAILED;
    }
    cblas_daxpy(n, 1.0, b, 1, r, 1);
 
    /* --------------------------------
     * Outer iteration
    -------------------------------- */
    while (it < itmax)
    {
        beta = cblas_dnrm2(n, r, 1);
        if (beta < tol) {
            *iter = it;
            return UCFD_STATUS_CONVERGED;
        }
        
        // Left-preconditioning
        psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, r);
 
        beta = cblas_dnrm2(n, r, 1);
        y[0] = beta;
 
        // V = r/beta
        #pragma omp parallel for
        for (i=0; i<n; i++) {
            V[i] = r[i]/beta;
        }
 
        /* --------------------------------
         * Inner iteration (Restart)
        -------------------------------- */
        for (j = 0; j < m; j++)
        {
            mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, op, descr, &V[j * n], 0.0, w);
            if (mklstat != SPARSE_STATUS_SUCCESS) {
                *iter = mklstat;
                return UCFD_MKL_FAILED;
            }
            psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, w);
 
            // Arnoldi iteration
            for (i = 0; i < (j + 1); i++) {
                tmp = cblas_ddot(n, w, 1, &V[i * n], 1);
                H[j + m * i] = tmp;
                cblas_daxpy(n, -tmp, &V[i * n], 1, w, 1);
            }
 
            tmp = cblas_dnrm2(n, w, 1);
            H[j + (j + 1) * m] = tmp;
            
            #pragma omp parallel for
            for (i=0; i<n; i++) {
                V[(j+1)*n+i] = w[i]/tmp;
            }
 
            // Givens Rotation
            for (i = 0; i < j; i++)
            {
                c = g[i * 2];     // g[i, 0]
                s = g[i * 2 + 1]; // g[i, 1]
                h1 = H[j + i * m];
                h2 = H[j + (i + 1) * m];
                H[j + i * m] = c * h1 - s * h2;
                H[j + (i + 1) * m] = s * h1 + c * h2;
            }
 
            h1 = H[j * (m + 1)];
            h2 = H[j + (j + 1) * m];
            rr = sqrt(h1 * h1 + h2 * h2);
            c = h1 / rr;
            s = -h2 / rr;
            H[j * (m + 1)] = c * h1 - s * h2;
            g[j * 2] = c;
            g[j * 2 + 1] = s;
 
            // Modify e1 vector
            y[j + 1] = y[j];
            y[j] *= c;
            y[j + 1] *= s;
        }
 
        // Back substitution
        cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit, m, H, m, y, 1);
 
        // Update
        cblas_dgemv(CblasRowMajor, CblasTrans, m, n, 1.0, V, n, y, 1, 1.0, x, 1);
 
        // Computes next iteration residual
        mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1.0, op, descr, x, 0.0, r);
        if (mklstat != SPARSE_STATUS_SUCCESS) {
            *iter = mklstat;
            return UCFD_MKL_FAILED;
        }
        cblas_daxpy(n, 1.0, b, 1, r, 1);
        
        ++it;
    }
    return UCFD_MAX_ITER;
}
 
ucfd_status_t step_gmres(sparse_matrix_t op, ucfd_precon_solve psolve, const struct matrix_descr descr,
                         int bn, int m, int *flag,
                         int *row_ptr, int *col_ind, int *diag_ind, double *precon_nnz_data,
                         double *x, double *b, double *H, double *V, double *g, double *y, double *w, double *r)
{
    int i, j;
    const int n = bn * BLOCK;
    double beta, tmp, c, s, h1, h2, rr;
    sparse_status_t mklstat;
    *flag = 0;
 
    // Apply preconditioner
    psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, r);
    
    beta = cblas_dnrm2(n, r, 1);
    y[0] = beta;
 
    // V = r/beta
    #pragma omp parallel for
    for (i=0; i<n; i++) {
        V[i] = r[i]/beta;
    }
 
    /* --------------------------------
     * Inner iteration (Restart)
    -------------------------------- */
    for (j = 0; j < m; j++)
    {
        mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, op, descr, &V[j * n], 0.0, w);
        if (mklstat != SPARSE_STATUS_SUCCESS) {
            *flag = mklstat;
            return UCFD_MKL_FAILED;
        }
 
        psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, w);
 
        // Arnoldi iteration
        for (i = 0; i < (j + 1); i++)
        {
            tmp = cblas_ddot(n, w, 1, &V[i * n], 1);
            H[j + m * i] = tmp;
            cblas_daxpy(n, -tmp, &V[i * n], 1, w, 1);
        }
 
        tmp = cblas_dnrm2(n, w, 1);
        H[j + (j + 1) * m] = tmp;
 
        #pragma omp parallel for
        for (i=0; i<n; i++) {
            V[(j+1)*n+i] = w[i]/tmp;
        }
 
        // Givens Rotation
        for (i = 0; i < j; i++)
        {
            c = g[i * 2];     // g[i, 0]
            s = g[i * 2 + 1]; // g[i, 1]
            h1 = H[j + i * m];
            h2 = H[j + (i + 1) * m];
            H[j + i * m] = c * h1 - s * h2;
            H[j + (i + 1) * m] = s * h1 + c * h2;
        }
 
        h1 = H[j * (m + 1)];
        h2 = H[j + (j + 1) * m];
        rr = sqrt(h1 * h1 + h2 * h2);
        c = h1 / rr;
        s = -h2 / rr;
        H[j * (m + 1)] = c * h1 - s * h2;
        g[j * 2] = c;
        g[j * 2 + 1] = s;
 
        // Modify e1 vector
        y[j + 1] = y[j];
        y[j] *= c;
        y[j + 1] *= s;
    }
 
    // Back substitution
    cblas_dtrsv(CblasRowMajor, CblasUpper, CblasNoTrans, CblasNonUnit, m, H, m, y, 1);
 
    // Update
    cblas_dgemv(CblasRowMajor, CblasTrans, m, n, 1.0, V, n, y, 1, 1.0, x, 1);
 
    // Computes next iteration residual
    mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1.0, op, descr, x, 0.0, r);
    if (mklstat != SPARSE_STATUS_SUCCESS) {
        *flag = mklstat;
        return UCFD_MKL_FAILED;
    }
    cblas_daxpy(n, 1.0, b, 1, r, 1);
 
    return UCFD_STATUS_SUCCESS;
}
 
ucfd_status_t serial_bicgstab(sparse_matrix_t op, ucfd_precon_type_t precon_type, int bn, int *iter, double tol,
                              int *row_ptr, int *col_ind, int *diag_ind, double *precon_nnz_data,
                              double *x, double *b, double *r, double *p, double *v, double *s, double *t)
{
    int it, itmax;
    const int n = bn * BLOCK;
    double rho, rhoprev, alpha, beta, omega, resid;
    double rv, ts, tt;
    ucfd_precon_solve psolve;
    sparse_status_t mklstat;
 
    // Sparse matrix description for system matrix A(op)
    struct matrix_descr descr;
    descr.type = SPARSE_MATRIX_TYPE_GENERAL;
    descr.mode = 0;
    descr.diag = 0;
 
    // Set preconditioner solver
    if (precon_type == BILU) psolve = bilu_psolve;
    else if (precon_type == LUSGS) psolve = lusgs_psolve;
    else psolve = none_psolve;
 
    itmax = *iter;
    it = 0;
 
    // Computes residual : r := b - A @ x
    mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1.0, op, descr, x, 0.0, r);
    if (mklstat != SPARSE_STATUS_SUCCESS) {
        *iter = mklstat;
        return UCFD_MKL_FAILED;
    }
    cblas_daxpy(n, 1.0, b, 1, r, 1);
 
    // Choose r\tilde as r
    /* r[:n] = r
    r[n:] = r\tilde */
    cblas_dcopy(n, r, 1, &r[n], 1);
 
    // Outer iteration
    while (it < itmax)
    {
        resid = cblas_dnrm2(n, r, 1);
        if (resid < tol)
            return UCFD_STATUS_CONVERGED;
 
        rho = cblas_ddot(n, r, 1, &r[n], 1);
 
        // Rho breakdown
        if (fabs(rho) < eps)
            return UCFD_STATUS_RHO_BREAKDOWN;
 
        if (it == 0)
            cblas_dcopy(n, r, 1, p, 1);
        else
        {
            beta = (rho / rhoprev) * (alpha / omega);
 
            // p = r + beta*(p - omega*v)
            cblas_daxpy(n, -omega, v, 1, p, 1);
            cblas_dscal(n, beta, p, 1);
            cblas_daxpy(n, 1.0, r, 1, p, 1);
        }
 
        // phat = inv(M) @ p
        cblas_dcopy(n, p, 1, &p[n], 1);
        psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, &p[n]);
 
        // v = A @ phat
        mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, op, descr, &p[0], 0.0, v);
        if (mklstat != SPARSE_STATUS_SUCCESS)
        {
            *iter = mklstat;
            return UCFD_MKL_FAILED;
        }
 
        rv = cblas_ddot(n, v, 1, &r[n], 1);
        alpha = rho / rv;
 
        // s = r - alpha*v
        // 1) r := r - alpha*v (r = s)
        cblas_daxpy(n, -alpha, v, 1, r, 1);
        // 2) s := r
        cblas_dcopy(n, r, 1, s, 1);
 
        // x := x + alpha*phat
        cblas_daxpy(n, alpha, &p[n], 1, x, 1);
        resid = cblas_dnrm2(n, s, 1);
        if (resid < tol)
            return UCFD_STATUS_CONVERGED;
 
        // shat = inv(M) @ s
        cblas_dcopy(n, s, 1, &s[n], 1);
        psolve(bn, row_ptr, col_ind, diag_ind, precon_nnz_data, &s[n]);
 
        mklstat = mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, op, descr, &s[0], 0.0, t);
        if (mklstat != SPARSE_STATUS_SUCCESS)
        {
            *iter = mklstat;
            return UCFD_MKL_FAILED;
        }
 
        ts = cblas_ddot(n, t, 1, s, 1);
        tt = cblas_ddot(n, t, 1, t, 1);
        omega = ts / tt;
 
        // Update solution
        cblas_daxpy(n, omega, &s[0], 1, x, 1);
 
        // r = s - omega*t
        cblas_daxpy(n, -omega, t, 1, r, 1);
 
        // Update rho
        rhoprev = rho;
 
        ++it;
    }
 
    return UCFD_MAX_ITER;
}