tscalapack.cpp

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

// STL
#include <iostream>

// MPI (must be before #define REAL)
#include "mpi.h"

using std::cout;
using std::endl;

extern "C" // {{{
{
    void   Cblacs_pinfo    (int* mypnum, int* nprocs);
    void   Cblacs_get      (int context, int request, int* value);
    int    Cblacs_gridinit (int* context, char * order, int np_row, int np_col);
    void   Cblacs_gridinfo (int context, int*  np_row, int* np_col, int*  my_row, int*  my_col);
    void   Cblacs_gridexit (int context);
    void   Cblacs_exit     (int error_code);
    int    numroc_         (int *n, int *nb, int *iproc, int *isrcproc, int *nprocs);
    void   descinit_       (int *desc, int *m, int *n, int *mb, int *nb, int *irsrc, int *icsrc, int *ictxt, int *lld, int *info);
    double pdlamch_        (int *ictxt , char *cmach);
    double pdlange_        (char *norm, int *m, int *n, double *A, int *ia, int *ja, int *desca, double *work);
    void   pdlacpy_        (char *uplo, int *m, int *n, double *a, int *ia, int *ja, int *desca, double *b, int *ib, int *jb, int *descb);
    void   pdgesv_         (int *n, int *nrhs, double *A, int *ia, int *ja, int *desca, int* ipiv, double *B, int *ib, int *jb, int *descb, int *info);
    void   pdgemm_         (char *TRANSA, char *TRANSB, int * M, int * N, int * K, double * ALPHA,
                            double * A, int * IA, int * JA, int * DESCA, double * B, int * IB, int * JB, int * DESCB,
                            double * BETA, double * C, int * IC, int * JC, int * DESCC );
    int    indxg2p_        (int *indxglob, int *nb, int *iproc, int *isrcproc, int *nprocs);
} // }}}

int main(int argc, char **argv)
{
    // Initialize
    MPI::Init(argc, argv);
    int rank = MPI::COMM_WORLD.Get_rank();
    int size = MPI::COMM_WORLD.Get_size();

    // Info
    cout << "Process #" << rank << " of " << size << " processes just started ...\n";
    cout << "" << endl;

    // Data
    int n     = 500; // nrows=ncols of a squared matrix
    int nrhs  = 1;   // number of Right-Hand-Side columns
    int nprow = 2;   // grid size
    int npcol = 2;   // grid size
    int nb    = 64;  // number of blocks

    // Check
    if (nb>n) nb=n;
    if (nprow*npcol>size)
    {
        if (rank==0)
        {
            cout << "ERROR: There is not enough processes (" << size << "<" << nprow*npcol << ") to make a p-by-q process grid" << endl;
            MPI::Finalize();
            return 1;
        }
    }

    // Initialize BLACS grid
    int iam,   nprocs, ictxt;
    int myrow, mycol;
    Cblacs_pinfo    (&iam, &nprocs);
    Cblacs_get      (-1, 0, &ictxt);
    Cblacs_gridinit (&ictxt, /*order*/"Col", nprow, npcol);
    Cblacs_gridinfo (ictxt, &nprow, &npcol, &myrow, &mycol);

    // Work only the process in the process grid
    if ((myrow>-1)&&(mycol>-1)&&(myrow<nprow)&&(mycol<npcol))
    {
        // Compute the size of the local matrices
        int izero = 0;
        int np    = numroc_(&n   , &nb, &myrow, &izero, &nprow);
        int nq    = numroc_(&n   , &nb, &mycol, &izero, &npcol);
        int nqrhs = numroc_(&nrhs, &nb, &mycol, &izero, &npcol);

        // Allocate and fill the matrices A and B
        double * A     = new double [np*nq];
        double * Acpy  = new double [np*nq];
        double * B     = new double [np*nqrhs];
        double * X     = new double [np*nqrhs];
        double * R     = new double [np*nqrhs];
        int    * ippiv = new int    [np+nb];

        // Generate random A matrix 
        int k=0;
        for (int i=0; i<np; i++)
        for (int j=0; j<nq; j++)
        {
            A[k] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
            k++;   
        }

        // Generate random B vector
        k=0;
        for (int i=0; i<np;    i++)
        for (int j=0; j<nqrhs; j++)
        {
            B[k] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
            k++;   
        }

        // Initialize the array descriptor for the matrix A and B
        int itemp = (np<1 ? 1 : np);
        int descA[9];
        int descB[9];
        int info;
        descinit_(descA, &n, &n   , &nb, &nb, &izero, &izero, &ictxt, &itemp, &info);
        descinit_(descB, &n, &nrhs, &nb, &nb, &izero, &izero, &ictxt, &itemp, &info);

        // Make a copy of A and the rhs for checking purposes
        int ione = 1;
        pdlacpy_("All", &n, &n   , A, &ione, &ione, descA, Acpy, &ione, &ione, descA);
        pdlacpy_("All", &n, &nrhs, B, &ione, &ione, descB, X   , &ione, &ione, descB);

        // Call ScaLAPACK PDGESV routine
        double MPIt1 = MPI::Wtime();
        pdgesv_(&n, &nrhs, A, &ione, &ione, descA, ippiv, X, &ione, &ione, descB, &info);
        double MPIt2 = MPI::Wtime();
        double MPIelapsed = MPIt2-MPIt1;

        // Froebenius norm = ||A * X  - B|| / ( ||X|| * ||A|| * eps * N )
        double mone=-1.0;
        double pone= 1.0;
        pdlacpy_("All", &n, &nrhs, B, &ione, &ione, descB, R, &ione, &ione, descB);
        pdgemm_ ("N", "N", &n, &nrhs, &n, &pone, Acpy, &ione, &ione, descA, X, &ione, &ione, descB, &mone, R, &ione, &ione, descB);
        double work   = 0.0;
        double eps    = pdlamch_(&ictxt, "Epsilon");
        double AnormF = pdlange_("F", &n, &n   , A, &ione, &ione, descA, &work);
        double XnormF = pdlange_("F", &n, &nrhs, X, &ione, &ione, descB, &work);
        double RnormF = pdlange_("F", &n, &nrhs, R, &ione, &ione, descB, &work);
        double residF = RnormF/(AnormF*XnormF*eps*static_cast<double>(n));

        // Information
        if (iam==0)
        {
            cout << endl;
            cout << "***********************************************" << endl;
            cout << "  Example of ScaLAPACK routine call: (PDGESV)  " << endl;
            cout << "***********************************************" << endl;
            cout << endl;
            cout << "  n      = " << n                      << endl;
            cout << "  nrhs   = " << nrhs                   << endl;
            cout << "  grid   = " << nprow << ", " << npcol << endl;
            cout << "  blocks = " << nb << "x" << nb        << endl;
            cout << "  np     = " << np                     << endl;
            cout << "  nq     = " << nq                     << endl;
            cout << "  nqrhs  = " << nqrhs                  << endl;
            cout << endl;
            cout << "MPI elapsed time during (pdgesv) = " << MPIelapsed << endl;
            cout << "Info returned by (pdgesv)        = " << info       << endl;
            cout << endl;
            cout << "Froebenius norm (residual) = " << residF << endl;
            cout << endl;
        }

        // Clean up
        delete [] A;
        delete [] Acpy;
        delete [] B;
        delete [] X;
        delete [] R;
        delete [] ippiv;

        // Grid exit
        Cblacs_gridexit(0);
    }

    // Information
    if(iam==0)
    {
        cout << "***********************************************\n";
        cout << "                      END                      \n";
        cout << "***********************************************\n";
    }

    // Finalize
    MPI::Finalize();

    return 0;
}

// vim:fdm=marker

---