# src/colarray/linalg.h

Go to the documentation of this file.
```00001 /***************************************************************************
00002                           linalg.h  -  description
00003                              -------------------
00004     begin                : Mon Dec 15 2003
00005     copyright            : (C) 2003 by Roman Geus
00006     email                : roman.geus@psi.ch
00007 ***************************************************************************/
00008
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017
00018 #ifndef LINALG_H
00019 #define LINALG_H
00020
00021 #include <cassert>
00022 #include <cmath>
00023 #include "vector.h"
00024 #include "matrix.h"
00025 #include "blas.h"
00026
00027 namespace linalg {
00028
00029     using namespace colarray;
00030
00031     template <typename T>
00032     T squared_sum(const Matrix<T>& A) {
00033         T sum = 0;
00034         for (int i = 0; i < A._m; i ++) {
00035             for (int j = 0; j < A._n; j ++) {
00036                 T val = A(i,j);
00037                 sum += val*val;
00038             }
00039         }
00040         return sum;
00041     }
00042
00043     long double norm_frobenius(const Matrix<long double>& A);
00044
00045     double norm_frobenius(const Matrix<double>& A);
00046
00047     float norm_frobenius(const Matrix<float>& A);
00048
00049     template <typename T>
00050     T determinant(const Matrix<T>& A) {
00051         assert(A._m == A._n);
00052
00053         T prod = 1;
00054         int sign = 1, info;
00055         Vector<int_t> pivots(A._m);
00056         Matrix<T> B(A._m, A._n);
00057
00058         B.inject(A);
00059         blas::getrf(B, pivots, info);
00060         assert(info >= 0);
00061         if (info == 0) {
00062             for (int i = 0; i < B._m; i ++) {
00063                 prod *= B(i,i);
00064                 if (pivots(i) != i+1)
00065                     sign = -sign;
00066             }
00067             return sign*prod;
00068         } else
00069             // U(info,info) == 0, singular matrix
00070             return static_cast<T>(0);
00071     }
00072
00073     template <typename T>
00074     void inverse(const Matrix<T>& A, Matrix<T>& B) {
00075         assert(A._m == A._n);
00076         assert(B._m == B._n);
00077         assert(A._m == B._m);
00078
00079         int n = A._m;
00080
00081         // B = idendity matrix
00082         B = 0;
00083         for (int i = 0; i < n; i ++)
00084             B(i,i) = 1;
00085
00086         // solve A*X = B, overwriting B with X
00087         int_t info;
00088         Vector<int_t> pivots(n);
00089         blas::gesv(A, B, pivots, info);
00090         assert(info == 0);
00091     }
00092
00095     template <typename T>
00096     void eigenvalues(const Matrix<T>& A, ColumnVector<T>& lambda) {
00097         assert(A._m == A._n);
00098         assert(A._m == lambda._n);
00099
00100         colarray::Matrix<T> B(A._m, A._n);
00101         B.inject(A);
00102         int_t info;
00103         blas::syev('n', 'u', B, lambda, info);
00104         assert(info == 0);
00105     }
00106
00107
00110     template <typename T>
00111       void eigenvalueDecomp(const Matrix<T>& A, Matrix<T>& Q, ColumnVector<T>& lambda) {
00112         assert(A._m == A._n);
00113         assert(Q._m == Q._n);
00114         assert(A._m == Q._m);
00115         assert(A._m == lambda._n);
00116
00117         Q.inject(A);
00118         int_t info;
00119         blas::syev('v', 'u', Q, lambda, info);
00120         assert(info == 0);
00121     }
00122
00123
00124 } // namespace linalg
00125
00126 #endif
```

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  1.4.7