home *** CD-ROM | disk | FTP | other *** search
- //===================================================================
- // linalg.h
- //
- // Version 1.1
- //
- // Written by:
- // Brent Worden
- // WordenWare
- // email: Brent@Worden.org
- //
- // Copyright (c) 1998-1999 WordenWare
- //
- // Created: August 28, 1998
- // Revised: April 17, 1999
- //===================================================================
-
- #ifndef _LINALG_H_
- #define _LINALG_H_
-
- #include "numerror.h"
- #include "matrix.hpp"
-
- NUM_BEGIN
-
- template<class T>
- bool choldcmp(Matrix<T>& a, Vector<T>& p)
- //-------------------------------------------------------------------
- // Given a positive-definite symmetric matrix a, this routine
- // constructs its Colesky decomposition, a = LL'. On input, only the
- // upper triangle of a need be given; it is not modified. The
- // Cholesky factor L is returned in the lower traingle of a, except
- // for its diagonal elements which are returned in p. The routine
- // returns true if the inversion was successful, otherwise it returns
- // false.
- //-------------------------------------------------------------------
- {
- int i, j, k;
- T sum;
- int n = a.rows();
-
- p[0] = sqrt(a[0][0]);
- for(i = 1; i < n; i++){
- a[i][0] /= p[0];
- }
- for(i = 1; i < n; ++i){
- sum = a[i][i];
- for(k = 0; k < i; ++k){
- sum -= a[i][k] * a[i][k];
- }
- if(sum <= 0.0){
- return false;
- }
- p[i] = sqrt(sum);
- for(j = n - 1; j > i; --j){
- sum = a[j][i];
- for(k = 0; k < i; ++k){
- sum -= a[j][k] * a[i][k];
- }
- a[j][i] = sum / p[i];
- }
- }
-
- return true;
- };
-
- template<class T>
- void cholsl(Matrix<T>& a, Vector<T>& p, Vector<T>& x, Vector<T>& b)
- //-------------------------------------------------------------------
- // Solves the set of n linear equations a.x = b, where a is a
- // positive-definite symmetric matrix. a and p are input as the
- // output of the routine choldcmp. Only the lower triangle of a is
- // accessed. b is input as the right-hand side vector. The solution
- // vector is returned in x[0..n-1]. a, and p are not modified and
- // can be left in place for successive calls with different
- // right-hand sides b. b is not modified unless you identify b as x
- // in the calling sequence, which is allowed.
- //-------------------------------------------------------------------
- {
- int i, k;
- T sum;
- int n = a.rows();
-
- for(i = 0; i < n; i++){
- for(sum = b[i], k = i - 1; k >= 0; --k){
- sum -= a[i][k] * x[k];
- }
- x[i] = sum / p[i];
- }
- for(i = n - 1; i >= 0; --i){
- for(sum = x[i], k = i + 1; k < n; ++k){
- sum -= a[k][i] * x[k];
- }
- x[i] = sum / p[i];
- }
- }
-
- template<class T>
- T determinant(Matrix<T>& a)
- //-------------------------------------------------------------------
- // Returns the determinant of the matrix a.
- //-------------------------------------------------------------------
- {
- Matrix<T> cp(a);
- T ret;
- int i;
- int n = cp.rows();
- Vector<int> indx(n);
- int sgn;
-
- if(ludcmp(cp, indx, sgn)){
- ret = T(sgn);
- for(i = 0; i < n; i++){
- ret *= cp[i][i];
- }
- } else {
- ret = 0.0;
- }
- return ret;
- };
-
- template<class T>
- bool inverse(Matrix<T>& a)
- //-------------------------------------------------------------------
- // Replaces the matrix a with its inverse. The routine returns true
- // if the inversion was successful, otherwise it returns false.
- //-------------------------------------------------------------------
- {
- int i, j;
- bool ret = false;
- int n = a.rows();
- Matrix<T> ai(n, n);
- Vector<int> indx(n);
- Vector<T> col(n);
- int d;
-
- if(ludcmp(a, indx, d)){
- for(j = 0; j < n; j++){
- for(i = 0; i < n; i++){
- col[i] = T(0);
- }
- col[j] = T(1);
- lubksb(a, indx, col);
- for(i = 0; i < n; i++){
- ai[i][j] = col[i];
- }
- }
- a = ai;
- ret = true;
- }
- return ret;
- };
-
- template<class T>
- void lubksb(Matrix<T>& a, Vector<int>& indx, Vector<T>& b)
- //-------------------------------------------------------------------
- // Given a matrix a and permutation vector indx returned from ludcmp,
- // this routines solves the set of linear equations a.x = b. On
- // input b holds the right-hand side vector. On output, it holds the
- // solution vector x. a and indx are not modified by this routine
- // and can be left in place for successive colls with different
- // right-hand sides b.
- //-------------------------------------------------------------------
- {
- int i, ii = -1, ip, j;
- T sum;
- int n = a.rows();
-
- for(i = 0; i < n; i++){
- ip = indx[i];
- sum = b[ip];
- b[ip] = b[i];
- if(ii > -1){
- for(j = ii; j <= i - 1; j++){
- sum -= a[i][j] * b[j];
- }
- } else if(sum){
- ii = i;
- }
- b[i] = sum;
- }
- for(i = n - 1; i >= 0; i--){
- sum = b[i];
- for(j = i + 1; j < n; j++){
- sum -= a[i][j] * b[j];
- }
- b[i] = sum / a[i][i];
- }
- };
-
- template<class T>
- bool ludcmp(Matrix<T>& a, Vector<int>& indx, int& d)
- //-------------------------------------------------------------------
- // See above comment. d is output as 1 or -1 depending on whether
- // the number of row interchanges was even or odd, respectively.
- // This routine is used in combination with lubksb to solve linear
- // equations or invert a matrix.
- //-------------------------------------------------------------------
- {
- int i, imax, j, k;
- T big, dum, sum, temp;
- int n = a.rows();
- Vector<T> vv(n);
- T zero(0);
- T one(1);
-
- d = 1;
- for(i = 0; i < n; i++){
- big = zero;
- for(j = 0; j < n; j++){
- if((temp = absoluteValue(a[i][j])) > big){
- big = temp;
- }
- }
- if(big == zero){
- throw Exception("num::ludcmp", "Singular Matrix");
- return false;
- }
- vv[i] = one / big;
- }
- for(j = 0; j < n; j++){
- for(i = 0; i < j; i++){
- sum = a[i][j];
- for(k = 0; k < i; k++){
- sum -= a[i][k] * a[k][j];
- }
- a[i][j] = sum;
- }
- big = zero;
- for(i = j; i < n; i++){
- sum = a[i][j];
- for(k = 0; k < j; k++){
- sum -= a[i][k] * a[k][j];
- }
- a[i][j] = sum;
- if((dum = vv[i] * absoluteValue(sum)) >= big){
- big = dum;
- imax = i;
- }
- }
- if(j != imax){
- for(k = 0; k < n; k++){
- dum = a[imax][k];
- a[imax][k] = a[j][k];
- a[j][k] = dum;
- }
- d = -d;
- vv[imax] = vv[j];
- }
- indx[j] = imax;
- if(a[j][j] == 0.0){
- a[j][j] = NUMERICS_FLOAT_MIN;
- }
- if(j != n - 1){
- dum = one / (a[j][j]);
- for(i = j + 1; i < n; i++){
- a[i][j] *= dum;
- }
- }
- }
- return true;
- };
-
- template<class T>
- bool linsolve(Matrix<T>& a, Vector<T>& b, Vector<T>& x)
- //-------------------------------------------------------------------
- // Solves the set of linear equations a.x = b. The routine returns
- // true if the solution vector is successfully found, otherwise it
- // returns false.
- //-------------------------------------------------------------------
- {
- Matrix<T> ac(a);
- Vector<T> bc(b);
- bool ret = false;
- int d;
- Vector<int> indx(a.rows());
-
- if(ludcmp(ac, indx, d)){
- lubksb(ac, indx, bc);
- x = bc;
- ret = true;
- }
- return ret;
- };
-
- NUM_END
-
- #endif
-
- //===================================================================
- // Revision History
- //
- // Version 1.0 - 08/28/1998 - New.
- // Version 1.1 - 04/10/1999 - Added Numerics namespace.
- // 04/17/1999 - Use matrix and vector classes.
- // 04/17/1999 - Removed dot product functions.
- //===================================================================
-