1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Dario Mangoni, Radu Serban
13 // =============================================================================
14
15 #include <mpi.h>
16 #include <bitset>
17
18 #include "chrono_mumps/ChMumpsEngine.h"
19
20 namespace chrono {
21
ChMumpsEngine()22 ChMumpsEngine::ChMumpsEngine() {
23 int myrank;
24 MPI_Init(nullptr, nullptr);
25 MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
26
27 /* Initialize a MUMPS instance. Use MPI_COMM_WORLD */
28 mumps_id.job = INIT;
29 mumps_id.par = 1;
30 mumps_id.sym = UNSYMMETRIC;
31 mumps_id.comm_fortran = -987654;
32
33 dmumps_c(&mumps_id);
34
35 /* Output messages */
36 mumps_id.ICNTL(1) = 6; // Error
37 mumps_id.ICNTL(2) = 0; // Diagnostic and warnings
38 mumps_id.ICNTL(3) = 0; // Global information
39 mumps_id.ICNTL(4) = 1; // Error, warning and diagnostic control
40
41 /* Matrix control */
42 mumps_id.ICNTL(5) = 0; // COO Matrix format selection
43 mumps_id.ICNTL(18) = 0; // Matrix centralized on the host
44 }
45
~ChMumpsEngine()46 ChMumpsEngine::~ChMumpsEngine() {
47 mumps_id.job = END;
48 dmumps_c(&mumps_id); // Terminate instance
49 MPI_Finalize();
50 }
51
SetProblem(const ChSparseMatrix & Z,ChVectorRef rhs)52 void ChMumpsEngine::SetProblem(const ChSparseMatrix& Z, ChVectorRef rhs) {
53 SetMatrix(Z);
54 SetRhsVector(rhs);
55 }
56
SetMatrix(const ChSparseMatrix & Z)57 void ChMumpsEngine::SetMatrix(const ChSparseMatrix& Z) {
58 // Convert to COO representation (1-indexed)
59 int dim = (int)Z.rows();
60 int nnz = (int)Z.nonZeros();
61
62 m_irn.resize(nnz);
63 m_jcn.resize(nnz);
64 m_a.resize(nnz);
65
66 int i = 0;
67 for (int k = 0; k < Z.outerSize(); ++k) {
68 for (ChSparseMatrix::InnerIterator it(Z, k); it; ++it) {
69 m_irn[i] = (int)it.row() + 1;
70 m_jcn[i] = (int)it.col() + 1;
71 m_a[i] = it.value();
72 i++;
73 }
74 }
75
76 mumps_id.n = dim;
77 mumps_id.nz = nnz;
78 mumps_id.irn = m_irn.data();
79 mumps_id.jcn = m_jcn.data();
80 mumps_id.a = m_a.data();
81 }
82
SetMatrixSymmetry(mumps_SYM mat_type)83 void ChMumpsEngine::SetMatrixSymmetry(mumps_SYM mat_type) {
84 mumps_id.sym = mat_type;
85 }
86
SetRhsVector(ChVectorRef b)87 void ChMumpsEngine::SetRhsVector(ChVectorRef b) {
88 mumps_id.rhs = b.data();
89 }
90
SetRhsVector(double * b)91 void ChMumpsEngine::SetRhsVector(double* b) {
92 mumps_id.rhs = b;
93 }
94
EnableNullPivotDetection(bool val,double threshold)95 void ChMumpsEngine::EnableNullPivotDetection(bool val, double threshold) {
96 mumps_id.ICNTL(24) = val; ///< activates null pivot detection
97 mumps_id.ICNTL(25) = 0; ///< tries to compute one of the many solutions of AX = B
98 mumps_id.CNTL(5) = 1e20; ///< fixation value
99 mumps_id.CNTL(3) = threshold; ///< pivot threshold
100 }
101
102
MumpsCall(mumps_JOB job_call)103 int ChMumpsEngine::MumpsCall(mumps_JOB job_call) {
104 /* Call the MUMPS package. */
105 mumps_id.job = job_call;
106 dmumps_c(&mumps_id);
107 return mumps_id.INFOG(1);
108 }
109
PrintINFOG()110 void ChMumpsEngine::PrintINFOG() {
111 if (mumps_id.INFOG(1) > 0) {
112 printf("WARN: INFOG(1)=");
113 typedef std::bitset<sizeof(int)> IntBits;
114 if (IntBits(mumps_id.INFOG(1)).test(0))
115 printf("+1: Row or column index out of range. Faulty entries: %d\n", mumps_id.INFOG(2));
116 if (IntBits(mumps_id.INFOG(1)).test(1))
117 printf("+2: Solution max-norm = 0\n");
118 if (IntBits(mumps_id.INFOG(1)).test(3))
119 printf("+8: More than %d iterative refinements are required\n", mumps_id.ICNTL(10));
120 return;
121 }
122
123 if (mumps_id.INFOG(1) == 0) {
124 printf("INFOG(1)=0: Mumps is successful!\n");
125 return;
126 }
127
128 printf("ERR: INFOG(1)=");
129 switch (mumps_id.INFOG(1)) {
130 case (0):
131 printf("0: Mumps is successful!\n");
132 break;
133 case (-1):
134 printf("-1: Error on processor %d\n", mumps_id.INFOG(2));
135 break;
136 case (-2):
137 printf("-2: Number of nonzeros out of range NZ=%d\n", mumps_id.INFOG(2));
138 break;
139 case (-3):
140 printf("-3: Mumps called with wrong JOB. JOB=%d\n", mumps_id.INFOG(2));
141 break;
142 case (-4):
143 printf("-4: Error in user-provided permutation array PERM_IN at position: %d\n", mumps_id.INFOG(2));
144 break;
145 case (-5):
146 printf("-5: Problem of real workspace allocation of size %d during analysis\n", mumps_id.INFOG(2));
147 break;
148 case (-6):
149 printf("-6: Matrix is singular in structure. Matrix rank: %d\n", mumps_id.INFOG(2));
150 break;
151 case (-7):
152 printf("-7: Problem of integer workspace allocation of size %d during analysis\n", mumps_id.INFOG(2));
153 break;
154 case (-10):
155 printf("-10: Matrix is numerically singular.\n");
156 break;
157 case (-16):
158 printf("-16: N is out of range. N=%d\n", mumps_id.INFOG(2));
159 break;
160 case (-21):
161 printf("-21: PAR=1 not allowed because only one processor is available.\n");
162 break;
163 case (-22):
164 printf("-22: Array pointers have problems. INFOG(2)=%d\n", mumps_id.INFOG(2));
165 break;
166 default:
167 printf("%d: See the user guide. INFOG(2)=%d\n", mumps_id.INFOG(1), mumps_id.INFOG(2));
168 }
169 }
170
171 } // namespace chrono
172