PLaSK library

Loading...
Searching...
No Matches
iterative_matrix.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK_COMMON_FEM_ITERATIVE_MATRIX_H
15#define PLASK_COMMON_FEM_ITERATIVE_MATRIX_H
16
17#include <algorithm>
18
19#include <nspcg/nspcg.hpp>
20
21#include "matrix.hpp"
22
23 namespace plask {
24
25 struct IterativeMatrixParams {
27 enum Accelelator {
28 ACCEL_CG,
29 ACCEL_SI,
30 ACCEL_SOR,
31 ACCEL_SRCG,
32 ACCEL_SRSI,
33 ACCEL_BASIC,
34 ACCEL_ME,
35 ACCEL_CGNR,
36 ACCEL_LSQR,
37 ACCEL_ODIR,
38 ACCEL_OMIN,
39 ACCEL_ORES,
40 ACCEL_IOM,
41 ACCEL_GMRES,
42 ACCEL_USYMLQ,
43 ACCEL_USYMQR,
44 ACCEL_LANDIR,
45 ACCEL_LANMIN,
46 ACCEL_LANRES,
47 ACCEL_CGCR,
48 ACCEL_BCGS
49 };
50 Accelelator accelerator = ACCEL_CG;
51
53 enum Preconditioner {
54 PRECOND_RICH,
55 PRECOND_JAC,
56 PRECOND_LJAC,
57 PRECOND_LJACX,
58 PRECOND_SOR,
59 PRECOND_SSOR,
60 PRECOND_IC,
61 PRECOND_MIC,
62 PRECOND_LSP,
63 PRECOND_NEU,
64 PRECOND_LSOR,
65 PRECOND_LSSOR,
66 PRECOND_LLSP,
67 PRECOND_LNEU,
68 PRECOND_BIC,
69 PRECOND_BICX,
70 PRECOND_MBIC,
71 PRECOND_MBICX
72 };
73 Preconditioner preconditioner = PRECOND_IC;
74
75 int maxit = 1000;
76 double maxerr = 1e-6;
77 int nfact = 10;
78
79 int ns1 = 5;
80 int ns2 = 100000;
81 int lvfill = 0;
82 int ltrunc = 0;
83 int ndeg = 1;
84 double omega = 1.0;
85
86 enum NoConvergenceBehavior { NO_CONVERGENCE_ERROR, NO_CONVERGENCE_WARNING, NO_CONVERGENCE_CONTINUE };
87 NoConvergenceBehavior no_convergence_behavior = NO_CONVERGENCE_WARNING;
89
90 // Output parameters
91 bool converged = true;
92 int iters = 0;
93 double err = 0;
94};
95
96 struct SparseMatrix : FemMatrix {
97 typedef void (*NspcgFunc)(...);
98
99 protected:
100 int* icords;
101
102 IterativeMatrixParams* params;
103
104 const int nstore, ndim, mdim;
105
106 int ifact = 1;
107 int nw = 0, inw = 0;
108 double* wksp = nullptr;
109 int* iwksp = nullptr;
110 int kblsz = -1, nbl2d = -1;
111
112 virtual NspcgFunc get_preconditioner() = 0;
113
114 virtual int get_maxnz() const { return mdim; }
115
116 public:
117 template <typename SolverT>
118 SparseMatrix(SolverT* solver, size_t rank, size_t size, size_t isiz)
119 : FemMatrix(solver, rank, size),
120 icords(aligned_malloc<int>(isiz)),
121 params(&solver->iter_params),
122 nstore(2),
123 ndim(rank),
124 mdim(isiz) {}
125
126 template <typename SolverT>
127 SparseMatrix(SolverT* solver, size_t rank, size_t size)
128 : FemMatrix(solver, rank, size),
129 icords(aligned_malloc<int>(2 * size)),
130 params(&solver->iter_params),
131 nstore(4),
132 ndim(size),
133 mdim(size) {}
134
135 ~SparseMatrix() {
136 aligned_free<int>(icords);
137 aligned_free<double>(wksp);
138 aligned_free<int>(iwksp);
139 }
140
141 void solverhs(DataVector<double>& B, DataVector<double>& X) override {
142 iparm_t iparm;
143 rparm_t rparm;
144 nspcg_dfault(iparm, rparm);
145
146 iparm.nstore = nstore;
147 iparm.itmax = params->maxit;
148 iparm.ipropa = 0;
149 iparm.ifact = (--ifact) ? 0 : 1;
150 if (ifact <= 0) ifact = params->nfact;
151 rparm.zeta = params->maxerr;
152
153 iparm.ns1 = params->ns1;
154 iparm.ns2 = params->ns2;
155 iparm.lvfill = params->lvfill;
156 iparm.ltrunc = params->ltrunc;
157 iparm.ndeg = params->ndeg;
158 rparm.omega = params->omega;
159
160 iparm.ns3 =
161 (params->accelerator == IterativeMatrixParams::ACCEL_LANMIN || params->accelerator == IterativeMatrixParams::ACCEL_CGCR)
162 ? 40
163 : 0;
164 iparm.kblsz = kblsz;
165 iparm.nbl2d = nbl2d;
166
167 solver->writelog(LOG_DETAIL, "Iterating linear system");
168
169#ifdef NDEBUG
170 iparm.level = -1;
171#else
172 iparm.level = 3;
173#endif
174
175 int maxnz = get_maxnz();
176 size_t default_nw = 3 * rank + 2 * params->maxit + rank * maxnz + std::max(kblsz, 1);
177 size_t default_inw = maxnz + std::max(2 * int(rank), maxnz * maxnz + maxnz);
178
179 if (nw < default_nw) {
180 nw = default_nw;
181 aligned_free<double>(wksp);
182 wksp = aligned_malloc<double>(nw);
183 iparm.ifact = 1; // we need to do factorization with new workspace
184 }
185
186 if (inw < default_inw) {
187 inw = default_inw;
188 aligned_free<int>(iwksp);
189 iwksp = aligned_malloc<int>(inw);
190 iparm.ifact = 1; // we need to do factorization with new workspace
191 }
192
193 // for (size_t r = 0; r < rank; ++r) {
194 // for (size_t c = 0; c < rank; ++c) {
195 // if (std::find(icords, A.icords+(ld+1), std::abs(int(r)-int(c))) == icords+(ld+1) )
196 // std::cout << " . ";
197 // else
198 // std::cout << str((*this))(r, c), "{:8.3f}") << " ";
199 // }
200 // std::cout << " " << str(B[r], "{:8.3f}") << std::endl;
201 // }
202
203 assert(B.size() == rank);
204
205 DataVector<double> U;
206 if (X.data() == nullptr || X.data() == B.data())
207 U.reset(B.size(), 1.);
208 else
209 U = X;
210 assert(U.size() == B.size());
211
212 int ier;
213
214 // TODO add choice of algorithms and predonditioners
215
216 NspcgFunc precond_func, accel_func;
217
218 if ((params->accelerator == IterativeMatrixParams::ACCEL_SOR) !=
219 (params->preconditioner == IterativeMatrixParams::PRECOND_SOR ||
220 params->preconditioner == IterativeMatrixParams::PRECOND_LSOR)) {
221 throw BadInput(solver->getId(), "SOR oraccelerator must be used with SOR or LSOR preconditioner");
222 }
223 if (params->accelerator == IterativeMatrixParams::ACCEL_SRCG &&
224 params->preconditioner != IterativeMatrixParams::PRECOND_SSOR &&
225 params->preconditioner != IterativeMatrixParams::PRECOND_LSSOR) {
226 throw BadInput(solver->getId(), "SRCG accelerator must be used with SSOR or LSSOR preconditioner");
227 }
228 if (params->accelerator == IterativeMatrixParams::ACCEL_SRSI &&
229 params->preconditioner != IterativeMatrixParams::PRECOND_SSOR &&
230 params->preconditioner != IterativeMatrixParams::PRECOND_LSSOR) {
231 throw BadInput(solver->getId(), "SRSI accelerator must be used with SSOR or LSSOR preconditioner");
232 }
233 // clang-format off
234 precond_func = this->get_preconditioner();
235
236 switch (params->accelerator) {
237 case IterativeMatrixParams::ACCEL_CG: accel_func = nspcg_cg; break;
238 case IterativeMatrixParams::ACCEL_SI: accel_func = nspcg_si; break;
239 case IterativeMatrixParams::ACCEL_SOR: accel_func = nspcg_sor; break;
240 case IterativeMatrixParams::ACCEL_SRCG: accel_func = nspcg_srcg; break;
241 case IterativeMatrixParams::ACCEL_SRSI: accel_func = nspcg_srsi; break;
242 case IterativeMatrixParams::ACCEL_BASIC: accel_func = nspcg_basic; break;
243 case IterativeMatrixParams::ACCEL_ME: accel_func = nspcg_me; break;
244 case IterativeMatrixParams::ACCEL_CGNR: accel_func = nspcg_cgnr; break;
245 case IterativeMatrixParams::ACCEL_LSQR: accel_func = nspcg_lsqr; break;
246 case IterativeMatrixParams::ACCEL_ODIR: accel_func = nspcg_odir; break;
247 case IterativeMatrixParams::ACCEL_OMIN: accel_func = nspcg_omin; break;
248 case IterativeMatrixParams::ACCEL_ORES: accel_func = nspcg_ores; break;
249 case IterativeMatrixParams::ACCEL_IOM: accel_func = nspcg_iom; break;
250 case IterativeMatrixParams::ACCEL_GMRES: accel_func = nspcg_gmres; break;
251 case IterativeMatrixParams::ACCEL_USYMLQ: accel_func = nspcg_usymlq; break;
252 case IterativeMatrixParams::ACCEL_USYMQR: accel_func = nspcg_usymqr; break;
253 case IterativeMatrixParams::ACCEL_LANDIR: accel_func = nspcg_landir; break;
254 case IterativeMatrixParams::ACCEL_LANMIN: accel_func = nspcg_lanmin; break;
255 case IterativeMatrixParams::ACCEL_LANRES: accel_func = nspcg_lanres; break;
256 case IterativeMatrixParams::ACCEL_CGCR: accel_func = nspcg_cgcr; break;
257 case IterativeMatrixParams::ACCEL_BCGS: accel_func = nspcg_bcgs; break;
258 };
259 // clang-format on
260
261 while (true) {
262 nspcg(precond_func, accel_func, ndim, mdim, rank, maxnz, data, icords, nullptr, nullptr, U.data(), nullptr, B.data(),
263 wksp, iwksp, nw, inw, iparm, rparm, ier);
264
265 // Increase workspace if needed
266 if (ier == -2 && nw) {
267 aligned_free<double>(wksp);
268 wksp = aligned_malloc<double>(nw);
269 iparm.ifact = 1; // we need to do factorization with new workspace
270 } else if (ier == -3 && inw) {
271 aligned_free<int>(iwksp);
272 iwksp = aligned_malloc<int>(inw);
273 iparm.ifact = 1; // we need to do factorization with new workspace
274 } else
275 break;
276 }
277
278 if (ier != 0) {
279 switch (ier) {
280 case -1: throw ComputationError(solver->getId(), "nonpositive matrix rank {}", rank);
281 case -2: throw ComputationError(solver->getId(), "insufficient real workspace ({} required)", nw);
282 case -3: throw ComputationError(solver->getId(), "insufficient integer workspace ({} required)", inw);
283 case -4: throw ComputationError(solver->getId(), "nonpositive diagonal element in stiffness matrix");
284 case -5: throw ComputationError(solver->getId(), "nonexistent diagonal element in stiffness matrix");
285 case -6: throw ComputationError(solver->getId(), "stiffness matrix A is not positive definite");
286 case -7: throw ComputationError(solver->getId(), "preconditioned matrix Q is not positive definite");
287 case -8: throw ComputationError(solver->getId(), "cannot permute stiffness matrix as requested");
288 case -9:
289 throw ComputationError(solver->getId(),
290 "Number of non-zero diagonals is not large enough to allow expansion of matrix");
291 case -10: throw ComputationError(solver->getId(), "inadmissible parameter encountered");
292 case -11: throw ComputationError(solver->getId(), "incorrect storage mode for block method");
293 case -12: throw ComputationError(solver->getId(), "zero pivot encountered in factorization");
294 case -13: throw ComputationError(solver->getId(), "breakdown in direction vector calculation");
295 case -14: throw ComputationError(solver->getId(), "breakdown in attempt to perform rotation");
296 case -15: throw ComputationError(solver->getId(), "breakdown in iterate calculation");
297 case -16: throw ComputationError(solver->getId(), "unimplemented combination of parameters");
298 case -18: throw ComputationError(solver->getId(), "unable to perform eigenvalue estimation");
299 case 1:
300 params->converged = false;
301 switch (params->no_convergence_behavior) {
302 case IterativeMatrixParams::NO_CONVERGENCE_ERROR:
303 throw ComputationError(solver->getId(), "failed to converge in {} iterations (error {})", iparm.itmax,
304 rparm.zeta);
305 case IterativeMatrixParams::NO_CONVERGENCE_WARNING:
306 solver->writelog(LOG_WARNING, "Failed to converge in {} iterations (error {})", iparm.itmax,
307 rparm.zeta);
308 break;
309 case IterativeMatrixParams::NO_CONVERGENCE_CONTINUE:
310 solver->writelog(LOG_DETAIL, "Did not converge yen in {} iterations (error {})", iparm.itmax,
311 rparm.zeta);
312 break;
313 }
314 break;
315 case 2: solver->writelog(LOG_WARNING, "`maxerr` was too small, reset to {}", 3.55e-12); break;
316 case 3:
317 solver->writelog(LOG_DEBUG,
318 "NSPGS: `zbrent` failed to converge in the maximum number of {} iterations (signifies "
319 "difficulty in eigenvalue estimation)",
320 std::max(params->maxit, 50));
321 break;
322 case 4:
323 solver->writelog(LOG_DEBUG,
324 "NSPGS: In `zbrent`, f (a) and f (b) have the same sign (signifies difficulty in "
325 "eigenvalue estimation)");
326 break;
327 case 5: solver->writelog(LOG_DEBUG, "NSPGS: Negative pivot encountered in factorization"); break;
328 }
329 }
330 if (ier != 1) {
331 solver->writelog(LOG_DETAIL, "Converged after {} iterations (error {})", iparm.itmax, rparm.zeta);
332 params->converged = true;
333 }
334
335 params->iters = iparm.itmax;
336 params->err = rparm.zeta;
337
338 if (X.data() != U.data()) X = U;
339 }
340
341 void mult(const DataVector<const double>& vector, DataVector<double>& result) override {
342 std::fill(result.begin(), result.end(), 0.);
343 addmult(vector, result);
344 }
345};
346
347 struct SparseBandMatrix : SparseMatrix {
354 template <typename SolverT>
355 SparseBandMatrix(SolverT* solver, size_t rank, size_t major) : SparseMatrix(solver, rank, 5 * rank, 5) {
356 icords[0] = 0;
357 icords[1] = 1;
358 icords[2] = major - 1;
359 icords[3] = major;
360 icords[4] = major + 1;
361 kblsz = major - 1;
362 nbl2d = major - 1;
363 }
364
372 template <typename SolverT>
373 SparseBandMatrix(SolverT* solver, size_t rank, size_t major, size_t minor) : SparseMatrix(solver, rank, 14 * rank, 14) {
374 icords[0] = 0;
375 icords[1] = 1;
376 icords[2] = minor - 1;
377 icords[3] = minor;
378 icords[4] = minor + 1;
379 icords[5] = major - minor - 1;
380 icords[6] = major - minor;
381 icords[7] = major - minor + 1;
382 icords[8] = major - 1;
383 icords[9] = major;
384 icords[10] = major + 1;
385 icords[11] = major + minor - 1;
386 icords[12] = major + minor;
387 icords[13] = major + minor + 1;
388 kblsz = minor - 1;
389 nbl2d = major - minor - 1;
390 }
391
398 template <typename SolverT>
399 SparseBandMatrix(SolverT* solver, size_t rank, std::initializer_list<int> bands)
400 : SparseMatrix(solver, rank, bands.size() * rank, bands.size()) {
401 size_t i = 0;
402 for (int band : bands) icords[i++] = band;
403 assert(icords[0] == 0);
404 }
405
412 double& operator()(size_t r, size_t c) override {
413 if (r == c) return data[r];
414 if (r < c) std::swap(r, c);
415 size_t i = std::find(icords, icords + mdim, r - c) - icords;
416 assert(i != mdim);
417 return data[c + rank * i];
418 }
419
420 void addmult(const DataVector<const double>& vector, DataVector<double>& result) override {
421 for (size_t r = 0; r < rank; ++r) {
422 result[r] += data[r] * vector[r];
423 }
424 for (size_t d = 1; d < mdim; ++d) {
425 size_t sd = rank * d;
426 for (size_t r = 0; r < rank; ++r) {
427 size_t c = r + icords[d];
428 if (c >= rank) break;
429 result[r] += data[r + sd] * vector[c];
430 result[c] += data[r + sd] * vector[r];
431 }
432 }
433 }
434
435 protected:
436 NspcgFunc get_preconditioner() override {
437 switch (params->preconditioner) {
438 case IterativeMatrixParams::PRECOND_RICH: return nspcg_rich2;
439 case IterativeMatrixParams::PRECOND_JAC: return nspcg_jac2;
440 case IterativeMatrixParams::PRECOND_LJAC: return nspcg_ljac2;
441 case IterativeMatrixParams::PRECOND_LJACX: return nspcg_ljacx2;
442 case IterativeMatrixParams::PRECOND_SOR: return nspcg_sor2;
443 case IterativeMatrixParams::PRECOND_SSOR: return nspcg_ssor2;
444 case IterativeMatrixParams::PRECOND_IC: return nspcg_ic2;
445 case IterativeMatrixParams::PRECOND_MIC: return nspcg_mic2;
446 case IterativeMatrixParams::PRECOND_LSP: return nspcg_lsp2;
447 case IterativeMatrixParams::PRECOND_NEU: return nspcg_neu2;
448 case IterativeMatrixParams::PRECOND_LSOR: return nspcg_lsor2;
449 case IterativeMatrixParams::PRECOND_LSSOR: return nspcg_lssor2;
450 case IterativeMatrixParams::PRECOND_LLSP: return nspcg_llsp2;
451 case IterativeMatrixParams::PRECOND_LNEU: return nspcg_lneu2;
452 case IterativeMatrixParams::PRECOND_BIC: return nspcg_bic2;
453 case IterativeMatrixParams::PRECOND_BICX: return nspcg_bicx2;
454 case IterativeMatrixParams::PRECOND_MBIC: return nspcg_mbic2;
455 case IterativeMatrixParams::PRECOND_MBICX: return nspcg_mbicx2;
456 };
457 assert(NULL);
458 return nullptr;
459 }
460
461 public:
462 void setBC(DataVector<double>& B, size_t r, double val) override {
463 data[r] = 1.;
464 B[r] = val;
465 // above diagonal
466 for (ptrdiff_t i = mdim - 1; i > 0; --i) {
467 ptrdiff_t c = r - icords[i];
468 if (c >= 0) {
469 ptrdiff_t ii = c + rank * i;
470 assert(ii < size);
471 B[c] -= data[ii] * val;
472 data[ii] = 0.;
473 }
474 }
475 // below diagonal
476 for (ptrdiff_t i = 1; i < mdim; ++i) {
477 ptrdiff_t c = r + icords[i];
478 if (c < rank) {
479 size_t ii = r + rank * i;
480 assert(ii < size);
481 B[c] -= data[ii] * val;
482 data[ii] = 0.;
483 }
484 }
485 }
486};
487
488 struct SparseFreeMatrix : SparseMatrix {
489 int inz;
490 int* const ir;
491 int* const ic;
492
499 template <typename SolverT>
500 SparseFreeMatrix(SolverT* solver, size_t rank, size_t maxnz)
501 : SparseMatrix(solver, rank, maxnz), inz(rank), ir(icords), ic(icords + maxnz) {
502 assert(maxnz >= rank);
503 for (size_t i = 0; i < rank; ++i) ir[i] = i + 1;
504 for (size_t i = 0; i < rank; ++i) ic[i] = i + 1;
505 if (params->preconditioner == IterativeMatrixParams::PRECOND_IC)
506 params->preconditioner = IterativeMatrixParams::PRECOND_JAC;
507 }
508
515 double& operator()(size_t r, size_t c) override {
516 if (r == c) return data[r];
517 if (r > c) std::swap(r, c);
518 assert(inz < size);
519 ir[inz] = r + 1;
520 ic[inz] = c + 1;
521 return data[inz++];
522 }
523
524 void clear() override {
525 std::fill_n(data, size, 0.);
526 inz = rank;
527 }
528
529 void addmult(const DataVector<const double>& vector, DataVector<double>& result) override {
530 for (size_t i = 0; i < rank; ++i) {
531 result[i] += data[i] * vector[i];
532 }
533 for (size_t i = rank; i < inz; ++i) {
534 size_t r = ir[i] - 1, c = ic[i] - 1;
535 result[r] += data[i] * vector[c];
536 result[c] += data[i] * vector[r];
537 }
538 }
539
540 protected:
541 NspcgFunc get_preconditioner() override {
542 switch (params->preconditioner) {
543 case IterativeMatrixParams::PRECOND_RICH: return nspcg_rich4;
544 case IterativeMatrixParams::PRECOND_JAC: return nspcg_jac4;
545 case IterativeMatrixParams::PRECOND_LSP: return nspcg_lsp4;
546 case IterativeMatrixParams::PRECOND_NEU: return nspcg_neu4;
547 default: throw NotImplemented("preconditioner not implemented for non-rectangular or masked mesh");
548 };
549 assert(NULL);
550 }
551
552 int get_maxnz() const override { return inz; }
553
554 public:
555 void setBC(DataVector<double>& B, size_t r, double val) override {
556 data[r] = 1e32;
557 B[r] = val * 1e32;
558 }
559};
560
561} // namespace plask
562
563#endif // PLASK_COMMON_FEM_ITERATIVE_MATRIX_H