Line data Source code
1 : /*
2 : * Copyright (c) 2014: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include "super_lu.h"
34 :
35 : // fix warning
36 : #undef TRUE
37 : #undef FALSE
38 : #include "slu_ddefs.h"
39 :
40 : namespace ug{
41 :
42 :
43 : class SuperLUImplementation : public IExternalSolverImplementation
44 : {
45 : bool m_bInited;
46 : std::vector<double> rhs, nzval;
47 : /* row permutations from partial pivoting */
48 : /* column permutation vector */
49 : std::vector<int> perm_r, perm_c, colind, rowptr;
50 :
51 : SuperMatrix SuperLU_A, SuperLU_L, SuperLU_U, SuperLU_B;
52 :
53 : SuperLUConfiguration &config;
54 :
55 : superlu_options_t options;
56 : trans_t trans;
57 :
58 :
59 : public:
60 : SuperLUImplementation(SuperLUConfiguration &_config)
61 0 : : m_bInited(false), config(_config),
62 0 : trans(NOTRANS)
63 : {
64 0 : SuperLU_A.Store = NULL;
65 : }
66 :
67 0 : ~SuperLUImplementation()
68 0 : {
69 0 : destroy();
70 0 : }
71 :
72 0 : void destroy()
73 : {
74 0 : if(m_bInited)
75 : {
76 : // DO NOT destroy SuperLU_A and AA since we supplied all pointers!!!
77 0 : Destroy_SuperMatrix_Store(&SuperLU_B);
78 : memset(&SuperLU_B, 0, sizeof(SuperMatrix));
79 0 : Destroy_SuperNode_Matrix(&SuperLU_L);
80 : memset(&SuperLU_L, 0, sizeof(SuperMatrix));
81 0 : Destroy_CompCol_Matrix(&SuperLU_U);
82 : memset(&SuperLU_U, 0, sizeof(SuperMatrix));
83 :
84 0 : if (SuperLU_A.Store)
85 0 : SUPERLU_FREE(SuperLU_A.Store);
86 :
87 0 : m_bInited = false;
88 : }
89 0 : }
90 :
91 : void get_options(superlu_options_t& opt)
92 : {
93 : //opt.PrintStat = config.bPrintStat ? YES : NO;
94 0 : opt.Equil = config.equil ? YES : NO;
95 0 : switch(config.colPerm)
96 : {
97 0 : case SuperLUConfiguration::CPT_NATURAL: opt.ColPerm = NATURAL; break;
98 0 : case SuperLUConfiguration::CPT_MMD_ATA: opt.ColPerm = MMD_ATA; break;
99 0 : case SuperLUConfiguration::CPT_MMD_AT_PLUS_A: opt.ColPerm = MMD_AT_PLUS_A; break;
100 0 : case SuperLUConfiguration::CPT_COLAMD: opt.ColPerm = COLAMD; break;
101 : }
102 : }
103 :
104 : void
105 0 : dgssvA(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r,
106 : SuperMatrix *L, SuperMatrix *U, SuperMatrix *B,
107 : SuperLUStat_t *stat, int *info )
108 : {
109 : //DNformat *Bstore;
110 :
111 : int lwork = 0;
112 :
113 : /* Set default values for some parameters */
114 : int panel_size; /* panel size */
115 : int relax; /* no of columns in a relaxed snodes */
116 : int permc_spec;
117 0 : trans = NOTRANS;
118 :
119 :
120 : /* Convert A to SLU_NC format when necessary. */
121 : SuperMatrix* AA = NULL; // A in SLU_NC format used by the factorization routine.
122 : bool bAANeedsFree = false;
123 0 : if ( A->Stype == SLU_NR ) {
124 0 : NRformat *Astore = (NRformat*) A->Store;
125 0 : AA = (SuperMatrix*) SUPERLU_MALLOC(sizeof(SuperMatrix));
126 0 : dCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz,
127 0 : (double*) Astore->nzval, Astore->colind, Astore->rowptr,
128 : SLU_NC, A->Dtype, A->Mtype);
129 0 : trans = TRANS;
130 : bAANeedsFree = true;
131 : } else {
132 0 : if ( A->Stype == SLU_NC ) AA = A;
133 : }
134 :
135 : /*
136 : * Get column permutation vector perm_c[], according to permc_spec:
137 : * permc_spec = NATURAL: natural ordering
138 : * permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
139 : * permc_spec = MMD_ATA: minimum degree on structure of A'*A
140 : * permc_spec = COLAMD: approximate minimum degree column ordering
141 : * permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
142 : */
143 0 : permc_spec = options->ColPerm;
144 0 : if ( permc_spec != MY_PERMC && options->Fact == DOFACT )
145 0 : get_perm_c(permc_spec, AA, perm_c);
146 :
147 :
148 0 : int* etree = intMalloc(A->ncol);
149 :
150 : SuperMatrix AC; /* Matrix postmultiplied by Pc */
151 0 : sp_preorder(options, AA, perm_c, etree, &AC);
152 :
153 0 : panel_size = sp_ienv(1);
154 0 : relax = sp_ienv(2);
155 :
156 : /* Compute the LU factorization of A. */
157 : #ifdef SUPERLU_6_EXPERIMENTAL
158 : GlobalLU_t glu;
159 0 : dgstrf(options, &AC, relax, panel_size, etree,
160 : NULL, lwork, perm_c, perm_r, L, U, &glu, stat, info);
161 : #else
162 : dgstrf(options, &AC, relax, panel_size, etree,
163 : NULL, lwork, perm_c, perm_r, L, U, stat, info);
164 : #endif
165 :
166 0 : Destroy_CompCol_Permuted(&AC);
167 0 : if (bAANeedsFree)
168 : {
169 0 : SUPERLU_FREE(AA->Store);
170 0 : SUPERLU_FREE(AA);
171 : }
172 0 : SUPERLU_FREE(etree);
173 0 : }
174 :
175 :
176 : void
177 : dgssvB(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r,
178 : SuperMatrix *L, SuperMatrix *U, SuperMatrix *B,
179 : SuperLUStat_t *stat, int *info )
180 : {
181 0 : dgstrs (trans, L, U, perm_c, perm_r, B, stat, info);
182 : }
183 :
184 :
185 :
186 0 : virtual bool init(const CPUAlgebra::matrix_type &A)
187 : {
188 0 : destroy();
189 : PROFILE_BEGIN_GROUP(SuperLU_Preprocess, "algebra SuperLU");
190 : typedef CPUAlgebra::matrix_type::const_row_iterator row_it;
191 : typedef CPUAlgebra::matrix_type::value_type value_type;
192 :
193 0 : if( A.num_rows() == 0 || A.num_cols() == 0) return true;
194 :
195 : size_t numRows, numCols;
196 0 : A.copy_crs(numRows, numCols, nzval, rowptr, colind);
197 0 : THROW_IF_NOT_EQUAL(numRows, numCols);
198 : size_t N = numRows;
199 : size_t nnz = nzval.size();
200 : //if(N > 40000 && nnz > 400000) { UG_LOG("SuperLU preprocess, N = " << N << ", nnz = " << nnz << "... "); }
201 :
202 0 : dCreate_CompRow_Matrix(&SuperLU_A, N, N, nnz, &nzval[0], &colind[0], &rowptr[0], SLU_NR, SLU_D, SLU_GE);
203 :
204 : //rhs = doubleMalloc(N);
205 0 : rhs.resize(N, 0.0);
206 0 : dCreate_Dense_Matrix(&SuperLU_B, N, 1, &rhs[0], N, SLU_DN, SLU_D, SLU_GE);
207 :
208 0 : perm_r.resize(N+1);
209 0 : perm_c.resize(N+1);
210 : /* Set the default input options. */
211 0 : set_default_options(&options);
212 :
213 : get_options(options);
214 :
215 : SuperLUStat_t stat;
216 0 : StatInit(&stat);
217 : int info;
218 :
219 0 : dgssvA(&options, &SuperLU_A, &perm_c[0], &perm_r[0], &SuperLU_L, &SuperLU_U, &SuperLU_B, &stat, &info);
220 :
221 : /*
222 : dgssv_check_info(info, N);
223 : if(config.bPrintStat)
224 : StatPrint(&stat);
225 : */
226 0 : StatFree(&stat);
227 : //if(N > 40000 && nnz > 400000) { UG_LOG("done.\n"); }
228 0 : m_bInited = true;
229 0 : return true;
230 : }
231 :
232 0 : void dgssv_check_info(int info, size_t N)
233 : {
234 0 : if(info > 0)
235 : {
236 0 : if(info < (int)N)
237 : {
238 0 : UG_THROW("ERROR in SuperLU: U(i,i) with i=" << info << "is exactly zero. The factorization has\
239 : been completed, but the factor U is exactly singular,\
240 : so the solution could not be computed.");
241 : }
242 : else
243 0 : { UG_THROW("ERROR in SuperLU: memory allocation failure");}
244 : }
245 0 : else if(info < 0)
246 : {
247 0 : UG_THROW("ERROR in SuperLU: info < 0 ???");
248 : }
249 0 : }
250 :
251 0 : virtual bool apply(CPUAlgebra::vector_type &c, const CPUAlgebra::vector_type &d)
252 : {
253 : PROFILE_BEGIN_GROUP(SuperLU_Apply, "algebra SuperLU");
254 : size_t N = c.size();
255 0 : if(N == 0) return true;
256 0 : double *b = (double*) ((DNformat*) SuperLU_B.Store)->nzval;
257 :
258 0 : for (size_t i = 0; i < N; ++i)
259 0 : b[i] = d[i];
260 :
261 : superlu_options_t options;
262 0 : set_default_options(&options);
263 0 : options.Fact = FACTORED;
264 : int info;
265 : SuperLUStat_t stat;
266 0 : StatInit(&stat);
267 0 : dgssvB(&options, &SuperLU_A, &perm_c[0], &perm_r[0], &SuperLU_L, &SuperLU_U, &SuperLU_B, &stat, &info);
268 0 : StatFree(&stat);
269 0 : dgssv_check_info(info, N);
270 :
271 0 : for (size_t i = 0; i < N; ++i)
272 0 : c[i] = b[i];
273 :
274 : return true;
275 : }
276 :
277 0 : virtual const char* name() const { return "SuperLU"; }
278 : };
279 :
280 :
281 0 : IExternalSolverImplementation *CreateSuperLUImplementation(SuperLUConfiguration &config)
282 : {
283 0 : return new SuperLUImplementation(config);
284 : }
285 :
286 : }
287 :
|
