1 #ifndef theplu_yat_utility_diagonal_matrix 2 #define theplu_yat_utility_diagonal_matrix 25 #include "BasicMatrix.h" 26 #include "BasicVector.h" 27 #include "BLAS_utility.h" 29 #include "MatrixExpression.h" 30 #include "VectorExpression.h" 33 #include "yat_assert.h" 76 size_t rows(
void)
const;
86 const double operator()(
size_t row,
size_t col)
const;
106 class DiagonalMatrixMatrix :
111 : lhs_(lhs), rhs_(rhs) {}
112 size_t rows(
void)
const {
return lhs_.rows(); }
113 size_t columns(
void)
const {
return rhs_.columns(); }
114 double operator()(
size_t i,
size_t j)
const 115 {
return lhs_(i, i) * rhs_(i, j); }
117 void calculate_matrix(gsl_matrix*& result)
const 119 detail::reallocate(result, rows(), columns());
120 for (
size_t i=0; i<rows(); ++i)
121 for (
size_t j=0; j<columns(); ++j)
122 gsl_matrix_set(result, i, j, (*
this)(i, j));
126 const DiagonalMatrix& lhs_;
127 const BasicMatrix<T>& rhs_;
132 class DiagonalMatrixVector
133 :
public VectorExpression<DiagonalMatrixVector<T> >
136 DiagonalMatrixVector(
const DiagonalMatrix& lhs,
const BasicVector<T>& rhs)
137 : m_(lhs), vec_(rhs), size_(lhs.rows())
141 DiagonalMatrixVector(
const BasicVector<T>& lhs,
const DiagonalMatrix& rhs)
142 : m_(rhs), vec_(lhs), size_(rhs.columns())
146 double operator()(
size_t i)
const 148 YAT_ASSERT(i < size());
149 return m_(i, i) * vec_(i);
153 size_t size(
void)
const 159 void calculate_gsl_vector_p(
void)
const 161 this->allocate_memory(size_);
162 for (
size_t i=0; i<size_; ++i)
163 gsl_vector_set(this->v_, i, (*
this)(i));
166 const DiagonalMatrix& m_;
167 const BasicVector<T>& vec_;
174 class MatrixDiagonalMatrix :
175 public MatrixExpression<MatrixDiagonalMatrix<T> >
178 MatrixDiagonalMatrix(
const BasicMatrix<T>& lhs,
const DiagonalMatrix& rhs)
179 : lhs_(lhs), rhs_(rhs) {}
180 size_t rows(
void)
const {
return lhs_.rows(); }
181 size_t columns(
void)
const {
return rhs_.columns(); }
182 double operator()(
size_t i,
size_t j)
const 183 {
return lhs_(i, j) * rhs_(j, j); }
185 void calculate_matrix(gsl_matrix*& result)
const 187 detail::reallocate(result, rows(), columns());
188 for (
size_t i=0; i<rows(); ++i)
189 for (
size_t j=0; j<columns(); ++j)
190 gsl_matrix_set(result, i, j, (*
this)(i, j));
194 const BasicMatrix<T>& lhs_;
195 const DiagonalMatrix& rhs_;
200 class VectorDiagonalMatrix :
201 public BasicVector<VectorDiagonalMatrix<T> >
204 VectorDiagonalMatrix(
const BasicVector<T>& rhs,
const DiagonalMatrix& lhs)
205 : lhs_(lhs), rhs_(rhs) {}
207 double operator()(
size_t i)
const 209 YAT_ASSERT(i < size());
210 return lhs_(i) * rhs_(i, i);
213 size_t size(
void)
const 215 return rhs_.columns();
219 const gsl_vector* gsl_vector_p(
void)
const 222 gsl_vector* vec = detail::create_gsl_vector(n);
223 for (
size_t i=0; i<n; ++i)
224 gsl_vector_set(vec, i, (*
this)(i));
229 const BasicVector<T>& lhs_;
230 const DiagonalMatrix& rhs_;
234 template<
class M1,
class M2,
class OP>
236 :
public MatrixExpression<DiagonalMatrix<M1, M2, OP> >
239 DiagonalMatrix(
const M1& lhs,
const M2& rhs)
240 : lhs_(lhs), rhs_(rhs)
242 YAT_ASSERT(lhs.rows() == rhs.rows());
243 YAT_ASSERT(lhs.columns() == rhs.columns());
245 double operator()(
size_t i,
size_t j)
const 247 return get(i, j, OP());
249 size_t rows(
void)
const {
return lhs_.rows(); }
250 size_t columns(
void)
const {
return rhs_.columns(); }
251 void calculate_matrix(gsl_matrix*& m)
const 253 detail::reallocate(m, rows(), columns());
254 for (
size_t i=0; i<rows(); ++i)
255 for (
size_t j=0; j<columns(); ++j)
256 gsl_matrix_set(m, i, j, (*
this)(i,j));
263 double get(
size_t i,
size_t j, Plus)
const 264 {
return lhs_(i, j) + rhs_(i, j); }
266 double get(
size_t i,
size_t j, Minus)
const 267 {
return lhs_(i, j) - rhs_(i, j); }
284 DiagonalMatrix
operator*(
const DiagonalMatrix& lhs,
285 const DiagonalMatrix& rhs);
295 expression::DiagonalMatrixMatrix<T>
298 YAT_ASSERT(lhs.columns() == rhs.
rows());
299 return expression::DiagonalMatrixMatrix<T>(lhs, rhs);
311 expression::MatrixDiagonalMatrix<T>
314 YAT_ASSERT(lhs.
columns() == rhs.rows());
315 return expression::MatrixDiagonalMatrix<T>(lhs, rhs);
338 expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>,
expression::Plus>
341 YAT_ASSERT(lhs.rows() == rhs.
rows());
342 YAT_ASSERT(lhs.columns() == rhs.
columns());
343 return expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>,
359 YAT_ASSERT(lhs.
rows() == rhs.rows());
360 YAT_ASSERT(lhs.
columns() == rhs.columns());
388 YAT_ASSERT(lhs.rows() == rhs.
rows());
389 YAT_ASSERT(lhs.columns() == rhs.
columns());
390 return expression::DiagonalMatrix<DiagonalMatrix, BasicMatrix<T>,
406 YAT_ASSERT(lhs.
rows() == rhs.rows());
407 YAT_ASSERT(lhs.
columns() == rhs.columns());
421 expression::DiagonalMatrixVector<T>
424 YAT_ASSERT(lhs.columns() == rhs.
size());
425 return expression::DiagonalMatrixVector<T>(lhs, rhs);
437 expression::DiagonalMatrixVector<T>
440 YAT_ASSERT(lhs.
size() == rhs.rows());
441 return expression::DiagonalMatrixVector<T>(lhs, rhs);
size_t size(void) const
Definition: BasicVector.h:71
expression::ScaledMatrix< T > operator*(const BasicMatrix< MatrixExpression< expression::ScaledMatrix< T > > > &A, double k)
Definition: BLAS_level3.h:162
expression::MatrixDiagonalMatrix< T > operator*(const BasicMatrix< T > &lhs, const DiagonalMatrix &rhs)
matrix matrix multiplication
Definition: DiagonalMatrix.h:312
The Department of Theoretical Physics namespace as we define it.
expression::DiagonalMatrixVector< T > operator*(const BasicVector< T > &lhs, const DiagonalMatrix &rhs)
vector matrix multiplication
Definition: DiagonalMatrix.h:438
size_t columns(void) const
Definition: BasicMatrix.h:67
expression::DiagonalMatrix< DiagonalMatrix, BasicMatrix< T >, expression::Minus > operator-(const DiagonalMatrix &lhs, const BasicMatrix< T > &rhs)
matrix matrix subtraction
Definition: DiagonalMatrix.h:386
expression::DiagonalMatrix< BasicMatrix< T >, DiagonalMatrix, expression::Minus > operator-(const BasicMatrix< T > &lhs, const DiagonalMatrix &rhs)
matrix matrix subtraction
Definition: DiagonalMatrix.h:404
This is the yat interface to GSL vector.
Definition: Vector.h:59
Definition: BLAS_utility.h:31
This is the yat interface to GSL vector.
Definition: VectorBase.h:55
const double operator()(size_t row, size_t col) const
size_t rows(void) const
Definition: BasicMatrix.h:61
An expression that can be converted to a Matrix.
Definition: MatrixExpression.h:46
Definition: BasicVector.h:48
Definition: BLAS_utility.h:30
expression::DiagonalMatrixMatrix< T > operator*(const DiagonalMatrix &lhs, const BasicMatrix< T > &rhs)
matrix matrix multiplication
Definition: DiagonalMatrix.h:296
Interface to GSL matrix.
Definition: Matrix.h:74
Definition: BasicMatrix.h:38
DiagonalMatrix(void)
Default constructor.
expression::DiagonalMatrixVector< T > operator*(const DiagonalMatrix &lhs, const BasicVector< T > &rhs)
matrix vector multiplication
Definition: DiagonalMatrix.h:422
expression::DiagonalMatrix< DiagonalMatrix, BasicMatrix< T >, expression::Plus > operator+(const DiagonalMatrix &lhs, const BasicMatrix< T > &rhs)
matrix matrix addition
Definition: DiagonalMatrix.h:339
Diagonal Matrix.
Definition: DiagonalMatrix.h:49
size_t columns(void) const
expression::ScaledMatrix< T > operator-(const BasicMatrix< MatrixExpression< expression::ScaledMatrix< T > > > &A)
Definition: BLAS_level3.h:194
expression::DiagonalMatrix< BasicMatrix< T >, DiagonalMatrix, expression::Plus > operator+(const BasicMatrix< T > &lhs, const DiagonalMatrix &rhs)
matrix matrix addition
Definition: DiagonalMatrix.h:357