examples/trilinos_mtx_demo.cpp

Go to the documentation of this file.

/***************************************************************************
                          main.cpp  -  description
                             -------------------
    begin                : Tue Dec  2 17:12:37 CET 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program 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.                                   *
 *                                                                         *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
//#include <config.h>
#endif

#include <sys/time.h>
#include <iostream>
#include <cstdlib>
#include <cassert>
#include <unistd.h>

#include "Epetra_config.h"

#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif

#include "AztecOO.h"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "matrixmarket/mtxreader0.h"
#include "matrixmarket/balancedmtxreader.h"

extern "C"
{
#include "pbe.h"
}

double second()
{
    struct timeval tv;

    int info = gettimeofday(&tv, NULL);
    assert(info == 0);
    return tv.tv_sec + tv.tv_usec / 1e6;
}

void usage(const Epetra_Comm& Comm, string err_msg)
{
    if (Comm.MyPID() == 0) {
        cerr << err_msg << endl << endl;
        cerr << "Usage: trilinos_mtx_demo [-r reader_type] mtx_file" << endl;
    }
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    exit(1);
}

int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm Comm;
#endif

#ifdef __GNUC__
    std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
#endif

    if (Comm.MyPID() == 0)
        cout << "running on " << Comm.NumProc() << " processor(s)..." << endl;

    extern char *optarg;
    extern int optind;

    int reader_type = 1; // default reader

    int c;
    while ((c = getopt(argc, argv, "r:")) != -1) {
        switch (c) {
        case 'r':
            reader_type = atoi(optarg);
            if (!(reader_type == 1 || reader_type == 2))
                usage(Comm, "Unsupported reader");
            break;
        case ':':
        case '?':
            usage(Comm, "Unsupported option");
            break;
        }
    }

    if (argc - optind != 1) {
        usage(Comm, "");
    }

    string filename(argv[optind]);

    if (access(filename.c_str(), R_OK)) {
        usage(Comm, "Unable to access matrix file.");
    }

    const char* app_name = NULL;
    BaseMtxReader *reader = NULL;
    switch (reader_type) {
    case 1:
        reader = new MtxReader0(filename, Comm);
        app_name = "trilinos_trivial";
        break;
    case 2:
        reader = new BalancedMtxReader(filename, Comm);
        app_name = "trilinos_balanced";
        break;
    }

    pbe_init(app_name, 4, 0, NULL);
    pbe_start(1, "whole_program");

    pbe_start(2, "matrix_input");
    double t = second();
    Epetra_CrsMatrix* pA = reader->read();
    Epetra_CrsMatrix& A = *pA;
    t = second() - t;
    pbe_stop(2);
    cout << "Matrix read and distribution time: " << t << " seconds." << endl;

    for (int p = 0; p < Comm.NumProc(); p ++) {
        if (Comm.MyPID() == p)
            cout << "processor " << p << ": " << A.NumMyNonzeros() << " non-zeros" << endl;
        Comm.Barrier();
    }

    // cout << A << endl;

    // right hand side
    Epetra_Vector b(A.RowMap());
    b.PutScalar(1.0);

    // initial guess
    Epetra_Vector x(A.RowMap());
    x.PutScalar(0.0);

    pbe_start(3, "solver_setup");
    t = second();
    // construct solver for A*x = b
    AztecOO::AztecOO solver(&A, &x, &b);
    solver.SetAztecOption(AZ_solver, AZ_bicgstab);
    solver.SetAztecOption(AZ_precond, AZ_Jacobi);
    solver.SetAztecOption(AZ_reorder, 1);
    //solver.SetAztecOption(AZ_subdomain_solve, AZ_ilu);
    t = second() - t;
    pbe_stop(3);
    cout << "Solver setup time: " << t << " seconds." << endl;

    t = second();

    // start iteration
    pbe_start(4, "iteration");
    int stat = solver.Iterate(2000, 1e-8);
    assert(!stat);
    pbe_stop(4);

    t = second() - t;
    cout << "Iteration time: " << t << " seconds." << endl;

    // calculate residual
    Epetra_Vector r(A.RowMap());
    A.Multiply(false, x, r);
    r.Update(1.0, b, -1.0);

    double norm_r;
    r.Norm2(&norm_r);

    if (Comm.MyPID() == 0)
        cout << "norm(b - A*x) = " << norm_r << endl;

    cout << "Solver performed " << solver.NumIters() << " iterations." << endl
         << "Norm of true residual = " << solver.TrueResidual() << endl;

    pbe_stop(1);
    pbe_dump();
    pbe_finalize(1);

    delete reader;

#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return EXIT_SUCCESS;
}

---