PLaSK library

Loading...
Searching...
No Matches
cholesky_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_CHOLESKY_MATRIX_H
15#define PLASK_COMMON_FEM_CHOLESKY_MATRIX_H
16
17#include <cstddef>
18
19#include "matrix.hpp"
20
21#define UPLO 'L'
22
23// BLAS routine to multiply matrix by vector
24#define dsbmv F77_GLOBAL(dsbmv, DSBMV)
25F77SUB dsbmv(const char& uplo,
26 const int& n,
27 const int& k,
28 const double& alpha,
29 const double* a,
30 const int& lda,
31 const double* x,
32 const int& incx,
33 const double& beta,
34 double* y,
35 const int& incy); // y = alpha*A*x + beta*y,
36
37// LAPACK routines to solve set of linear equations
38#define dpbtrf F77_GLOBAL(dpbtrf, DPBTRF)
39F77SUB dpbtrf(const char& uplo, const int& n, const int& kd, double* ab, const int& ldab, int& info);
40
41#define dpbtrs F77_GLOBAL(dpbtrs, DPBTRS)
42F77SUB dpbtrs(const char& uplo,
43 const int& n,
44 const int& kd,
45 const int& nrhs,
46 double* ab,
47 const int& ldab,
48 double* b,
49 const int& ldb,
50 int& info);
51
52 namespace plask {
53
58 struct DpbMatrix : BandMatrix {
64 DpbMatrix(const Solver* solver, size_t rank, size_t band)
65 : BandMatrix(solver, rank, band, ((band + 1 + (15 / sizeof(double))) & ~size_t(15 / sizeof(double))) - 1) {}
66
67 size_t index(size_t r, size_t c) {
68 assert(r < rank && c < rank);
69 if (r < c) {
70 assert(c - r <= kd);
71// if UPLO = 'U', AB(kd+i-j,j) = A(i,j) for max(0,j-kd)<=i<=j;
72// if UPLO = 'L', AB(i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
73#if UPLO == 'U'
74 return ld * c + r + kd;
75#else
76 return ld * r + c;
77#endif
78 } else {
79 assert(r - c <= kd);
80#if UPLO == 'U'
81 return ld * r + c + kd;
82#else
83 return ld * c + r;
84#endif
85 }
86 }
87
88 double& operator()(size_t r, size_t c) override { return data[index(r, c)]; }
89
90 void factorize() override {
91 int info = 0;
92
93 solver->writelog(LOG_DETAIL, "Factorizing system");
94
95 dpbtrf(UPLO, int(rank), int(kd), data, int(ld + 1), info);
96 if (info < 0)
97 throw CriticalException("{0}: Argument {1} of `dpbtrf` has illegal value", solver->getId(), -info);
98 else if (info > 0)
99 throw ComputationError(solver->getId(), "leading minor of order {0} of the stiffness matrix is not positive-definite",
100 info);
101 }
102
103 void solverhs(DataVector<double>& B, DataVector<double>& X) override {
104 solver->writelog(LOG_DETAIL, "Solving matrix system");
105
106 int info = 0;
107 dpbtrs(UPLO, int(rank), int(kd), 1, data, int(ld + 1), B.data(), int(B.size()), info);
108 if (info < 0) throw CriticalException("{0}: Argument {1} of `dpbtrs` has illegal value", solver->getId(), -info);
109
110 std::swap(B, X);
111 }
112
113 void mult(const DataVector<const double>& vector, DataVector<double>& result) override {
114 dsbmv(UPLO, int(rank), int(kd), 1.0, data, int(ld) + 1, vector.data(), 1, 0.0, result.data(), 1);
115 }
116
117 void addmult(const DataVector<const double>& vector, DataVector<double>& result) override {
118 dsbmv(UPLO, int(rank), int(kd), 1.0, data, int(ld) + 1, vector.data(), 1, 1.0, result.data(), 1);
119 }
120};
121
122} // namespace plask
123
124#endif // PLASK_COMMON_FEM_CHOLESKY_MATRIX_H