src/matrix_matrix.cpp

Go to the documentation of this file.

//
// C++ Implementation: matrix_matrix
//
// Description: 
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <algorithm>
#include <vector>
#include <set>
#include "rlog/rlog.h"
#include "Epetra_Comm.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_FECrsMatrix.h"
#include "Epetra_Export.h"
#include "Epetra_Import.h"
#include "Epetra_Map.h"
#include "matrix_matrix.h"

#define USE_MM_COLMAP 1

using namespace std;

static void bext_map_1(const Epetra_CrsMatrix& A, std::vector<int>& col_indices) {
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        int info;
        int num_entries_A;
        double *values_A;
        int *indices_A;

        // extract non-zeros of row i
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // loop over non-zero entries of row i of A
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);

            // add j_global to list of indices
            std::vector<int>::iterator found_it;
            found_it = std::find(col_indices.begin(),
                                 col_indices.end(),
                                 j_global);
            if (found_it == col_indices.end()) {
                col_indices.push_back(j_global);
            }
        }

    } // for i_local
}

static void bext_map_2(const Epetra_CrsMatrix& A, std::vector<int>& col_indices) {
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        int info;
        int num_entries_A;
        double *values_A;
        int *indices_A;

        // extract non-zeros of row i
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // loop over non-zero entries of row i of A
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            col_indices.push_back(j_global);
        }

    } // for i_local
    std::sort(col_indices.begin(), col_indices.end());
    std::vector<int>::iterator erase_it = std::unique(col_indices.begin(), col_indices.end());
    col_indices.erase(erase_it, col_indices.end());
}

static void bext_map_3(const Epetra_CrsMatrix& A, std::vector<int>& col_indices) {
    std::set<int> set_B;
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        int info;
        int num_entries_A;
        double *values_A;
        int *indices_A;

        // extract non-zeros of row i
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // loop over non-zero entries of row i of A
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            set_B.insert(j_global);
        }

    } // for i_local

    col_indices.reserve(set_B.size());
    for (std::set<int>::iterator it = set_B.begin();
        it != set_B.end();
        ++ it)
    {
        col_indices.push_back(*it);
    }
}

int poor_man_matrix_matrix(const Epetra_CrsMatrix& A, const Epetra_CrsMatrix& B, Epetra_CrsMatrix& C) {
    int info;

    if (A.Comm().NumProc() == 1) {
        // call optimised code for np == 1
        return poor_man_matrix_matrix_local(A, B, C);
    }

    // Calculate new map for B
    // Note: the indices need not to be sorted
#ifdef USE_MM_COLMAP
    Epetra_Map map_B(A.ColMap());
#else
    // Compute indices of rows of B needed
    std::vector<int> indices_B;
    bext_map_3(A, indices_B);
    rDebug("%d: poor_man_matrix_matrix: local rows: %d -> %d",
           A.Comm().MyPID(), B.NumMyRows(), indices_B.size());
    rDebug("%d: poor_man_matrix_matrix: column map size: %d",
           A.Comm().MyPID(), A.ColMap().NumMyElements());
    Epetra_Map map_B(-1, indices_B.size(), &indices_B[0], 0, A.Comm());
#endif

    // Redistribute B using new map
    Epetra_Import plan_B(map_B, B.RowMap());
    Epetra_CrsMatrix Bext(Copy, map_B, 0);
    info = Bext.Import(B, plan_B, Insert);
    rAssert(info == 0);
    info = Bext.FillComplete(B.DomainMap(), B.RangeMap());
    rAssert(info == 0);

    // loop over rows of A
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        // extract non-zeros of row i
        int num_entries_A, num_entries_B;
        double *values_A, *values_B;
        int *indices_A, *indices_B;
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // do local computation
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            rAssert(Bext.MyGRID(j_global));
            info = Bext.ExtractMyRowView(Bext.LRID(j_global), num_entries_B, values_B, indices_B);
            rAssert(info == 0);
            for (int l = 0; l < num_entries_B; ++ l) {
                double value = values_A[k] * values_B[l];
                int index = Bext.GCID(indices_B[l]);
                info = C.SumIntoGlobalValues(A.GRID(i_local), 1, &value, &index);
                if (info > 0)
                    info = C.InsertGlobalValues(A.GRID(i_local), 1, &value, &index);
                rAssert(info >= 0);
            } // for l
        } // for k

    } // for i_local
    info = C.FillComplete(B.OperatorDomainMap(), A.OperatorRangeMap());
    rAssert(info >= 0);
    return 0;
}

int poor_man_matrix_matrix_local(const Epetra_CrsMatrix& A, const Epetra_CrsMatrix& B, Epetra_CrsMatrix& C) {
    int info;

    // loop over rows of A
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        // extract non-zeros of row i
        int num_entries_A, num_entries_B;
        double *values_A, *values_B;
        int *indices_A, *indices_B;
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // do local computation
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            if (B.MyGRID(j_global)) {
                // row j_global of M is local
                info = B.ExtractMyRowView(B.LRID(j_global), num_entries_B, values_B, indices_B);
                rAssert(info == 0);
                for (int l = 0; l < num_entries_B; ++ l) {
                    double value = values_A[k] * values_B[l];
                    int index = B.GCID(indices_B[l]);
                    info = C.SumIntoGlobalValues(A.GRID(i_local), 1, &value, &index);
                    if (info > 0)
                        info = C.InsertGlobalValues(A.GRID(i_local), 1, &value, &index);
                    rAssert(info >= 0);
                } // for l
            }
        } // for k

    } // for i_local
    info = C.FillComplete(B.OperatorDomainMap(), A.OperatorRangeMap());
    rAssert(info >= 0);
    return 0;
}

int poor_man_matrix_matrix_graph(const Epetra_CrsGraph& A, const Epetra_CrsGraph& B, Epetra_CrsGraph& C) {
    int max_nnz_B = B.MaxNumIndices();
    vector<int> ind_B(max_nnz_B);
    int *indices_B = &ind_B[0];

    // loop over rows of A
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        // extract non-zeros of row i
        int num_entries_A, num_entries_B;
        int *indices_A;
        int info;
        info = A.ExtractMyRowView(i_local, num_entries_A, indices_A);
        rAssert(info == 0);

        // do local computation
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            if (B.MyGRID(j_global)) {
                // row j_global of M is local
                info = B.ExtractMyRowCopy(B.LRID(j_global), max_nnz_B, num_entries_B, indices_B);
                rAssert(info == 0);
                // transform to global column indices
                for (int l = 0; l < num_entries_B; ++ l)
                    indices_B[l] = B.GCID(indices_B[l]);
                // insert into C
                info = C.InsertGlobalIndices(A.GRID(i_local), num_entries_B, indices_B);
                rAssert(info >= 0);
            }
        } // for k
    } // for i_local
    return 0;
}

int poor_man_matrix_matrix_local2(const Epetra_CrsMatrix& A, const Epetra_CrsMatrix& B, Epetra_CrsMatrix& C) {
    int info;

    int max_nnz_B = B.MaxNumEntries();
    vector<int> ind_B(max_nnz_B);
    vector<double> val_B(max_nnz_B);
    double *values_B = &val_B[0];
    int *indices_B = &ind_B[0];

    // FIXME: A and C must have the same RowMap() or this routine fails!
    // FIXME: B and C must have the same ColMap() or this routine fails!

    // loop over rows of A
    for (int i_local = 0; i_local < A.NumMyRows(); ++ i_local) {
        // extract non-zeros of row i
        int num_entries_A, num_entries_B;
        double *values_A;
        int *indices_A;
        info = A.ExtractMyRowView(i_local, num_entries_A, values_A, indices_A);
        rAssert(info == 0);

        // do local computation
        for (int k = 0; k < num_entries_A; ++ k) {
            int j_global = A.GCID(indices_A[k]);
            if (B.MyGRID(j_global)) {
                // extract row j of B (having local indices)
                info = B.ExtractMyRowCopy(B.LRID(j_global), max_nnz_B, num_entries_B, values_B, indices_B);
                rAssert(info == 0);
                // scale extracted row by A[i,j]
                for (int l = 0; l < num_entries_B; ++ l)
                    values_B[l] *= values_A[k];
                // add scaled row into row i of C
                info = C.SumIntoMyValues(i_local, num_entries_B, values_B, indices_B);
                if (info != 0)
                    rDebug("info=%d", info);
                rAssert(info == 0);
            }
        } // for k

    } // for i_local
    info = C.FillComplete(B.OperatorDomainMap(), A.OperatorRangeMap());
    rAssert(info >= 0);
    return 0;
}

int transpose_matrix(const Epetra_CrsMatrix& A, Epetra_FECrsMatrix& T) {
    int info;
    
    for (int lrow = 0; lrow < A.NumMyRows(); lrow ++) {
        int grow = A.RowMap().GID(lrow);
        int nof_entries;
        int* col;
        double* val;
        info = A.ExtractMyRowView(lrow, nof_entries, val, col);
        rAssert(info == 0);
    
        for (int k = 0; k < nof_entries; k ++) {
            int gcol = A.ColMap().GID(col[k]);
            info = T.InsertGlobalValues(1, &gcol, 1, &grow, val+k);
            rAssert(info >= 0);
        }
    }
    
    info = T.GlobalAssemble(false);
    rAssert(info == 0);
    info = T.FillComplete(A.OperatorRangeMap(), A.OperatorDomainMap());
    rAssert(info == 0);
    return 0;   
}

---