Quantum Construct (qC++): src/qclinearalgebra.cpp Source File

00001 /***************************************************************************
00002  *   Quantum Construct (qC++)                                                                  *
00003  *   The Quantum Physics Computational Library                                         *
00004  *   Copyright (C) 2005 by Shekhar S. Chandra                                      *
00005  *   Shekhar.Chandra@sci.monash,edu.au                                     *
00006  *                                                                         *
00007  *   This program is free software; you can redistribute it and/or modify  *
00008  *   it under the terms of the GNU General Public License as published by  *
00009  *   the Free Software Foundation; either version 2 of the License, or     *
00010  *   (at your option) any later version.                                   *
00011  *                                                                         *
00012  *   This program is distributed in the hope that it will be useful,       *
00013  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00014  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00015  *   GNU General Public License for more details.                          *
00016  *                                                                         *
00017  *   You should have received a copy of the GNU General Public License     *
00018  *   along with this program; if not, write to the                         *
00019  *   Free Software Foundation, Inc.,                                       *
00020  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00021  ***************************************************************************/
00022 
00023 #include "qclinearalgebra.h"
00024 
00025 QCLinearAlgebra::QCLinearAlgebra()
00026 {
00027     //ctor
00028 }
00029 
00030 QCLinearAlgebra::~QCLinearAlgebra()
00031 {
00032     //dtor
00033 }
00034 
00035 void QCLinearAlgebra::solveTridiagonal(Array<parameterType,1>& a, Array<parameterType,1>& b, Array<parameterType,1>& c, Array<parameterType,1>& unknowns, Array<parameterType,1>& knowns)
00036 {
00037     Array<parameterType,1> tmp;
00038         int n = a.extent(firstDim), startElement = 1;
00039         parameterType bet;
00040 
00041         tmp.resize(n);
00042         if(b(startElement) == 0.0)
00043         {
00044                 cerr << "Tridiagonal Error - First Element of Diagonal is Zero" << endl;
00045                 exit(1);
00046         }
00047         //Decomposition and Forward Substitution
00048         unknowns(startElement) = knowns(startElement)/(bet = b(startElement));
00049         for(int j = startElement+1; j < n; j ++)
00050         {
00051                 tmp(j) = c(j-1)/bet;
00052                 bet = b(j) - a(j)*tmp(j);
00053                 if(bet == 0.0)
00054                 {
00055                         cerr << "Tridiagonal Error - Zero pivot found (check if Matrix A is Diagonally Dominant)" << endl;
00056                         exit(1);
00057                 }
00058                 else
00059                         unknowns(j) = (knowns(j) - a(j)*unknowns(j-1))/bet;
00060         }
00061         //Backsubstitution
00062         for(int j = n-2; j >= startElement; j --)
00063                 unknowns(j) = unknowns(j) - tmp(j+1)*unknowns(j+1);
00064 }
00065 
00066 void QCLinearAlgebra::solveTridiagonal(Array<complexType,1>& a, Array<complexType,1>& b, Array<complexType,1>& c, Array<complexType,1>& unknowns, Array<complexType,1>& knowns)
00067 {
00068         Array<complexType,1> tmp;
00069         int n = a.extent(firstDim), startElement = 1;
00070         complexType bet;
00071 
00072         tmp.resize(n);
00073         if(b(startElement) == complexZero)
00074         {
00075                 cerr << "Tridiagonal Error - First Element of Diagonal is Zero" << endl;
00076                 exit(1);
00077         }
00078         //Decomposition and Forward Substitution
00079         unknowns(startElement) = knowns(startElement)/(bet = b(startElement));
00080         for(int j = startElement+1; j < n; j ++)
00081         {
00082                 tmp(j) = c(j-1)/bet;
00083                 bet = b(j) - a(j)*tmp(j);
00084                 if(bet == complexZero)
00085                 {
00086                         cerr << "Tridiagonal Error - Zero pivot found (check if Matrix A is Diagonally Dominant)" << endl;
00087                         exit(1);
00088                 }
00089                 else
00090                         unknowns(j) = (knowns(j) - a(j)*unknowns(j-1))/bet;
00091         }
00092         //Backsubstitution
00093         for(int j = n-2; j >= startElement; j --)
00094                 unknowns(j) = unknowns(j) - tmp(j+1)*unknowns(j+1);
00095 }
00096 
00097 void QCLinearAlgebra::solveTridiagonal2(Array<complexType,1>& a, Array<complexType,1>& b, Array<complexType,1>& c, Array<complexType,1>& unknowns, Array<complexType,1>& knowns)
00098 {
00099         Array<complexType,1> tmp;
00100         int n = a.extent(firstDim), startElement = 1;
00101         tmp.resize(n);
00102         for(int j = startElement+1; j < n; j++)
00103         {
00104                 tmp(j) = a(j)/b(j-1);
00105                 b(j) = b(j) - tmp(j)*c(j-1);
00106                 knowns(j) = knowns(j) - tmp(j)*knowns(j-1);
00107         }
00108         unknowns(n-1) = knowns(n-1)/b(n-1);
00109         for(int j = n-2; j >= startElement; j--)
00110                 unknowns(j) = (knowns(j) - c(j)*unknowns(j+1))/b(j);
00111 }