/* Copyright (C) 2006-2009 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/ The program reads the matrix and RHS in the Matrix Market format then run MUMPS. This software is a copyrighted work licensed under the terms, described in the file "FREE_LICENSE". */ #include #include #include "matrices.h" #include "dmumps_c.h" // MUMPS does not have matrix vector product void mulMatrixByVec(int n, int symmetric, vector &a, vector &irn, vector &jcn, double *in, double *out) { for (int i = 0; i < n; ++i) out[i] = 0.; for (int k = 0; k < a.size(); k++) { out[irn[k] - 1] += a[k]*in[jcn[k] - 1]; if (symmetric && irn[k] != jcn[k]) out[jcn[k] - 1] += a[k]*in[irn[k] - 1]; } } int main(int argc, char **argv) { if (argc < 2) { cout << "run_mumps matrix rhs" << endl; return 1; } try { // 1) Reading matrix SparseMatrix mat; ifstream in(argv[1]); if (!in) { cout << "cannot open file " << argv[1] << endl; return 2; } Timing t0; mat.read(in); t0.write("Reading matrix"); if (mat.NRows() != mat.NCols()) { cout << "Matrix is not squared " << endl; return 3; } // 2) converting matrix to the MUMPS representation and setting up MUMPS // 2a) setting up MUMPS parameters and initialize it DMUMPS_STRUC_C id; id.job = -1; id.par = 1; if (mat.IsSymmetric) // if you know that you matrix is positive definite, id.sym should be 1, then it is faster // id.sym = 1; id.sym = 2; else id.sym = 0; id.comm_fortran = 1; Timing t1; dmumps_c(&id); t1.write("Initialization"); //streams id.icntl[0] = 6; id.icntl[1] = 6; id.icntl[2] = 6; id.icntl[3] = 3; //ordering metis (5), or pord (4), or AMD (0), AMF (2), QAMD (6) id.icntl[6] = 5; vector irn; vector jcn; vector a; // 2b) converting SparseMatrix to the MUMPS representation Timing r1; id.n = mat.NRows(); id.nz = mat.nnz; irn.resize(id.nz); jcn.resize(id.nz); a.resize(id.nz); int k = 0; for (int i = 0; i < id.n; ++i) for (map::iterator j = mat[i].begin(); j != mat[i].end(); ++j) { irn[k] = i + 1; jcn[k] = (*j).first + 1; a[k] = (*j).second; ++k; } id.irn = &*irn.begin(); id.jcn = &*jcn.begin(); id.a = &*a.begin(); mat.clearMemory(); r1.write("Converting matrix"); Timing t2; // 3) LU factorization id.job = 1; dmumps_c(&id); t2.write("Analysis"); Timing t3; id.job = 2; dmumps_c(&id); t3.write("Factorization"); if (argc < 3) return 0; // 4) Reading RHS and setting up MUMPS ifstream in2(argv[2]); if (!in2) { cout << "cannot open file " << argv[2] << endl; return 4; } Matrix rhs; rhs.read(in2); if (id.n != rhs.NRows()) { cout << "RHS is not compatible with the matrix" << endl; return 4; } Matrix rhs_copy(rhs); id.rhs = &*rhs.begin(); id.nrhs = rhs.NCols(); id.lrhs = id.n; Timing t4; // 5) Back substitution id.job = 3; dmumps_c(&id); t4.write("Back substitution"); // 6) compute the norm and write the solution out //solution in rhs (it is overwritten), the correct rhs is in rhs_copy vector out(id.n); for (int i = 0; i < rhs.NCols(); ++i) { mulMatrixByVec(id.n, id.sym, a, irn, jcn, rhs.column(i), &*out.begin()); double sum = 0.; for (int j = 0; j < id.n; ++j) sum += pow(rhs_copy(j, i) - out[j], 2); cout << "norm Ax - b for RHS " << i + 1 << " is " << sqrt(sum) << endl; } cout << "Solution is in " << argv[2] << ".solution" << endl; ofstream outf((string(argv[2]) + ".solution").c_str()); rhs.write(outf); // 7) clean up Timing t5; id.job = -2; dmumps_c(&id); t5.write("Cleaning MUMPS"); } catch (BaseError &t) { cout << "Error is encountered:" << endl; cout << t << endl; return 5; } catch (...) { cout << "Not known error is encountered." << endl; return 6; } return 0; }