examples/ml_maxwell_test.cpp

Go to the documentation of this file.

/***************************************************************************
    ml_maxwell_test.cpp  -  test program for Trilinos ML in Maxwell mode
                             -------------------
    begin                : Tue Apr 27 2004
    copyright            : (C) 2004 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.                                   *
 *                                                                         *
 ***************************************************************************/

// #include "Epetra_config.h"
#include <ml_config.h>

#ifdef HAVE_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif

// for the ML Maxwell preconditioner...
#if !defined(HAVE_ML_EPETRA) || !defined(HAVE_ML_TEUCHOS)
#error Trilinos needs to be configured with --enable-epetra --enable-teuchos
#endif

#include <ml_include.h>
#include <ml_epetra_preconditioner.h>
#include <Teuchos_ParameterList.hpp>

#include <AztecOO.h>
#include <Epetra_Operator.h>
#include <Epetra_MultiVector.h>
#include <Epetra_Vector.h>
#include <Epetra_LinearProblem.h>
#include "Epetra_Export.h"

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <unistd.h>

#include "matrixmarket/mtxreader0.h"
#include "matrixmarket/mtxreader.h"

void finalize() {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
}

int permute(const Epetra_Map& newHMap,
            const Epetra_Map& newKMap,
            Epetra_CrsMatrix*& K, 
            Epetra_CrsMatrix*& H, 
            Epetra_CrsMatrix*& Y) {

    Epetra_Export HExport(H->RangeMap(), newHMap);
    Epetra_CrsMatrix *Hnew = new Epetra_CrsMatrix(Copy, newHMap, 0);
    Hnew->Export(*H, HExport, Insert);
    Hnew->FillComplete(newHMap, newHMap);
    delete H;
    H = Hnew;

    Epetra_Export KExport(K->RangeMap(), newKMap);
    Epetra_CrsMatrix *Knew = new Epetra_CrsMatrix(Copy, newKMap, 0);
    Knew->Export(*K, KExport, Insert);
    Knew->FillComplete(newKMap, newKMap);
    delete K;
    K = Knew;

    Epetra_Export YExport(Y->RangeMap(), newKMap);
    Epetra_CrsMatrix *Ynew = new Epetra_CrsMatrix(Copy, newKMap, 0);
    Ynew->Export(*Y, YExport, Insert);
    Ynew->FillComplete(H->DomainMap(), newKMap);
    delete Y;
    Y = Ynew;

    return(0);
}

int permute_ml(Epetra_CrsMatrix*& K, Epetra_CrsMatrix*& H, Epetra_CrsMatrix*& Y) {
    const Epetra_Comm& Comm(K->Comm());
    const Epetra_Map& newHMap = H->RangeMap();
    const Epetra_Map& newMMap = K->RangeMap();

    int num_my_edges = newMMap.NumMyElements();
    int num_my_nodes = newHMap.NumMyElements();
    
    Epetra_Export YExport(Y->RangeMap(), newMMap);
    Epetra_CrsMatrix *Ynew = new Epetra_CrsMatrix(Copy, newMMap, 0);
    Ynew->Export(*Y, YExport, Insert);
    Ynew->FillComplete(newHMap, newMMap);
    delete Y;
    Y = Ynew;

    int num_my_elements = num_my_edges;
    int my_pid = Comm.MyPID();
    int num_proc = Comm.NumProc();

    int* local_size = new int[num_proc];
    Comm.GatherAll(&num_my_elements, local_size, 1);

    int i = 0;
    int j = 0;
    int maxs = 0;
    for (i = 0; i < num_proc; i ++)
        if (local_size[i] > maxs) 
            maxs = local_size[i]; 

    int* global_element = new int[maxs * num_proc];
    int* next_location = new int[maxs * num_proc];
    int offset = my_pid * maxs;

    int* col = new int[2];
    int* proc = new int[2];
    int* loc = new int[2];
    int num_entries = 0;
    int *indices = 0;
    double *values = 0;
    bool notfound = false;

    Epetra_Util util;
    util.SetSeed(763206768);

    //
    //  my replacement ("location") of the Epetra_Map.RemoteIDList method
    //
    int* local_s = new int[num_proc];
    local_s[my_pid] = num_my_nodes;
    Comm.GatherAll(local_s + my_pid, local_s, 1);
    int maxhs = 0;
    for (i = 0; i < num_proc; i ++)
        if (local_s[i] > maxhs) 
            maxhs = local_s[i];
    int* global_e = new int[maxhs * num_proc];
    int* location = new int[newHMap.NumGlobalElements()];
    int offs = my_pid * maxhs;
    for (i = 0; i < num_my_nodes; i ++) {
        global_e[offs + i] = newHMap.GID(i);
    }
    Comm.GatherAll(global_e + offs, global_e, maxhs);
    for (i = 0; i < num_proc; i++)
        for (j = i * maxhs; j < i * maxhs + local_s[i]; j ++) 
            location[global_e[j]] = i; 
    delete[] local_s;
    delete[] global_e;
    //
    // end 
    //

    for (i = 0; i < num_my_elements; i ++) {
        global_element[offset + i] = newMMap.GID(i);
        assert( Y->ExtractMyRowView(i, num_entries, values, indices) >= 0 );
        j = 0;
        notfound = true;
        while ((j < num_entries) && notfound) {
            col[j] = Y->GCID(indices[j]);

            //            newHMap.RemoteIDList(1, col + j, proc + j, loc + j);
            proc[j] = location[col[j]];

            if (proc[j] == my_pid) 
                notfound = false;
            else
                j ++;
        }
        if (notfound && (num_entries != 0)) {
            if (num_entries == 1)
                j = 0;
            else 
                j = util.RandomInt() % 2;
            next_location[offset + i] = proc[j];
        }
        else {
            next_location[offset + i] = my_pid; 
        }
    }
    for (i = offset + num_my_elements; i < offset + maxs; i ++)
        next_location[i] = -1; 
 
    Comm.GatherAll(next_location + offset, next_location, maxs);
    Comm.GatherAll(global_element + offset, global_element, maxs);    

    num_my_elements = 0;
    for (i = 0; i < maxs*num_proc; i ++) {
        if ( next_location[i] == my_pid )
            global_element[num_my_elements++] = global_element[i];
    }

    Epetra_Map MMap(newMMap.NumGlobalElements(), num_my_elements, global_element, 0, Comm);

    delete[] local_size;
    delete[] global_element;
    delete[] next_location;
    delete[] col;
    delete[] proc;
    delete[] loc;
    delete[] location;

    permute(newHMap, MMap, K, H, Y);

    return(0);
}


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;

    // register atexit function
    atexit(finalize);

    //
    // read and distribute MatrixMarket matrices
    //
    BaseMtxReader *reader = 0;

    // edge matrix
    string filename("K.mtx");
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access edge matrix file.\n";
        exit(1);
    }
    reader = new MtxReader0(filename, Comm);
    Epetra_CrsMatrix* pK = reader->read();
    delete reader;

    // node matrix
    filename = "H.mtx";
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access node matrix file.\n";
        exit(1);
    }
    reader = new MtxReader0(filename, Comm);
    Epetra_CrsMatrix* pH = reader->read();
    delete reader;

    // transfer matrix (with compatible maps)
    const Epetra_Map& edge_map = pK->DomainMap();
    const Epetra_Map& node_map = pH->RangeMap();
    filename = "Y.mtx";
    if (access(filename.c_str(), R_OK)) {
        cerr << "Unable to access transfer matrix file.\n";
        exit(1);
    }
    reader = new MtxReader(filename, edge_map, node_map, edge_map);
    Epetra_CrsMatrix* pY = reader->read();
    delete reader;

    // Permute matrices to make sure that edges are assigned to a processor, which also
    // holds at least one of its nodes.
    permute_ml(pK, pH, pY);

    //
    // Construct ML Maxwell preconditioner
    //
    Epetra_Operator *prec = 0;
    Teuchos::ParameterList MLList;

#ifndef SKIP_ML
    // set defaults for maxwell type problems
    ML_Epetra::SetDefaults("maxwell", MLList);
    MLList.set("aggregation: type", "Uncoupled-MIS");
        
    // Hack to bypass a typo in ML_Epetra::SetDefaults
    MLList.set("coarse: type", "Amesos-KLU");

    if (Comm.MyPID() == 0)
        MLList.print(cout, 8);
    
    // construct ML preconditioner
    prec = new ML_Epetra::MultiLevelPreconditioner(*pK,
                                                   *pY,
                                                   *pH,
                                                   MLList,
                                                   true);
#endif
    //
    // Solve linear system K * x = b
    //

    // Right hand side
    Epetra_Vector b(pK->OperatorRangeMap());
    b.PutScalar(1.0);

    // Initial guess
    Epetra_Vector x(pK->OperatorDomainMap());
    x.PutScalar(0.0);

    // Setup iterative solver
    Epetra_LinearProblem problem(pK, &x, &b);
    AztecOO solver(problem);
    solver.SetAztecOption(AZ_solver, AZ_bicgstab);

#ifndef SKIP_ML
    solver.SetPrecOperator(prec);
#else
    solver.SetAztecOption(AZ_precond, AZ_none);
#endif

    // Solve
    solver.Iterate(1000, 1.0E-8);

    // Compute Residual
    Epetra_Vector bcomp(pK->OperatorRangeMap());
    pK->Apply(x, bcomp);
    Epetra_Vector resid(pK->OperatorRangeMap());
    resid.Update(1.0, b, -1.0, bcomp, 0.0);
    double residual;
    resid.Norm2(&residual);
    if (Comm.MyPID() == 0)
        cout << "Residual    = " << residual << endl;

    // Free memory
    delete prec;
    delete pK;
    delete pH;
    delete pY;

    return EXIT_SUCCESS;
}

---