PLearn 0.1
BandedSolvers.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // BandedSolvers.h
00004 //
00005 // Copyright (C) 2004 Nicolas Chapados 
00006 // 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 /* *******************************************************      
00036  * $Id: BandedSolvers.h 9436 2008-09-04 18:48:55Z nouiz $ 
00037  ******************************************************* */
00038 
00039 // Authors: Nicolas Chapados
00040 
00044 #ifndef BandedSolvers_INC
00045 #define BandedSolvers_INC
00046 
00047 // Put includes here
00048 #include "TVec.h"
00049 
00050 namespace PLearn {
00051 using namespace std;
00052 
00076 template <class T>
00077 void PentadiagonalSolveInPlace(const TVec<T>& y, const TVec<T>& a,
00078                                const TVec<T>& b, const TVec<T>& c)
00079 {
00080     int n = a.size();
00081     int m = y.size();
00082     int o = b.size();
00083     int p = c.size();
00084 
00085     // Limitation of this routine for now
00086     if (m < 2)
00087         PLERROR("PentadiagonalSolve: vectors must have length at least two");
00088   
00089     if (n != m || n != o || n != p)
00090         PLERROR("PentadiagonalSolve: vector dimensions don't agree; they must "
00091                 "all have the same length.");
00092 
00093     c[m-1] = 0;
00094     c[m-2] = 0;
00095     b[m-1] = 0;
00096 
00097     T h1=0;
00098     T h2=0;
00099     T h3=0;
00100     T h4=0;
00101     T h5=0;
00102     T hh1=0;
00103     T hh2=0;
00104     T hh3=0;
00105     T hh5=0;
00106     T z=0;
00107     T hb=0;
00108     T hc=0;
00109 
00110     for (int i=0 ; i<m ; ++i) {
00111         z=a[i]-h4*h1-hh5*hh2;
00112         hb=b[i];
00113         hh1=h1;
00114         h1=(hb-h4*h2)/z;
00115         b[i]=h1;
00116         hc=c[i];
00117         hh2=h2;
00118         h2=hc/z;
00119         c[i]=h2;
00120         a[i]=(y[i]-hh3*hh5-h3*h4)/z;
00121         hh3=h3;
00122         h3=a[i];
00123         h4=hb-h5*hh1;
00124         hh5=h5;
00125         h5=hc;
00126     }
00127 
00128     h2=0;
00129     h1=a[m-1];
00130     y[m-1]=h1;
00131 
00132     for (int i=m-1 ; i>=0 ; --i) {
00133         y[i]=a[i]-b[i]*h1-c[i]*h2;
00134         h2=h1;
00135         h1=y[i];
00136     }
00137 }
00138 
00139 } // end of namespace PLearn
00140 
00141 #endif
00142 
00143 
00144 /*
00145   Local Variables:
00146   mode:c++
00147   c-basic-offset:4
00148   c-file-style:"stroustrup"
00149   c-file-offsets:((innamespace . 0)(inline-open . 0))
00150   indent-tabs-mode:nil
00151   fill-column:79
00152   End:
00153 */
00154 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines