PLearn 0.1
|
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 :