src/colarray/linalg.h

Go to the documentation of this file.

/***************************************************************************
                          linalg.h  -  description
                             -------------------
    begin                : Mon Dec 15 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef LINALG_H
#define LINALG_H

#include <cassert>
#include <cmath>
#include "vector.h"
#include "matrix.h"
#include "blas.h"

namespace linalg {

    using namespace colarray;

    template <typename T>
    T squared_sum(const Matrix<T>& A) {
        T sum = 0;
        for (int i = 0; i < A._m; i ++) {
            for (int j = 0; j < A._n; j ++) {
                T val = A(i,j);
                sum += val*val;
            }
        }
        return sum;
    }

    long double norm_frobenius(const Matrix<long double>& A);

    double norm_frobenius(const Matrix<double>& A);
  
    float norm_frobenius(const Matrix<float>& A);

    template <typename T>
    T determinant(const Matrix<T>& A) {
        assert(A._m == A._n);

        T prod = 1;
        int sign = 1, info;
        Vector<int_t> pivots(A._m);
        Matrix<T> B(A._m, A._n);

        B.inject(A);
        blas::getrf(B, pivots, info);
        assert(info >= 0);
        if (info == 0) {
            for (int i = 0; i < B._m; i ++) {
                prod *= B(i,i);
                if (pivots(i) != i+1)
                    sign = -sign;
            }
            return sign*prod;
        } else
            // U(info,info) == 0, singular matrix
            return static_cast<T>(0);
    }

    template <typename T>
    void inverse(const Matrix<T>& A, Matrix<T>& B) {
        assert(A._m == A._n);
        assert(B._m == B._n);
        assert(A._m == B._m);

        int n = A._m;

        // B = idendity matrix
        B = 0;
        for (int i = 0; i < n; i ++)
            B(i,i) = 1;

        // solve A*X = B, overwriting B with X
        int_t info;
        Vector<int_t> pivots(n);
        blas::gesv(A, B, pivots, info);
        assert(info == 0);
    }

    template <typename T>
    void eigenvalues(const Matrix<T>& A, ColumnVector<T>& lambda) {
        assert(A._m == A._n);
        assert(A._m == lambda._n);
    
        colarray::Matrix<T> B(A._m, A._n);
        B.inject(A);
        int_t info;
        blas::syev('n', 'u', B, lambda, info);
        assert(info == 0);
    }


    template <typename T>
      void eigenvalueDecomp(const Matrix<T>& A, Matrix<T>& Q, ColumnVector<T>& lambda) {
        assert(A._m == A._n);
        assert(Q._m == Q._n);
        assert(A._m == Q._m);
        assert(A._m == lambda._n);
    
        Q.inject(A);
        int_t info;
        blas::syev('v', 'u', Q, lambda, info);
        assert(info == 0);
    }


} // namespace linalg

#endif

---