! BoostConv -- Copyright 2010 Paolo Luchini ! http://CPLcode.net/Applications/Numerical/BoostConv/ !( Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. !) USE rbmat BoostConvWS=STRUCTURE[INTEGER CONSTANT N,lo,hi; ! BoostConv workspace keeps internal state between calls ! N,lo and hi will be specified at run time; see infocpl "variable size" INTEGER rot,rot2 ! indices used to circulate basis vectors ARRAY(0..N-1,0..N-1) OF REAL dd ! symmetric square matrix holding the scalar products of base vectors ARRAY(0..N-1,lo..hi) OF REAL in,out ! in and out base vectors of the boostconv algorithm ] POINTER TO BoostConvWS defaultBoostConvWS=NULL SUBROUTINE BoostConv[REAL r^(∗); OPTIONAL BoostConvWS ws^^=defaultBoostConvWS; INTEGER length=10] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Expedites the convergence of a pre-existing numerical iteration by ! constructing a probable guess of the forthcoming next residual as a ! linear combination of past residuals. ! ! Calling parameters: ! ! ws: pointer to internal workspace of type BoostConvWS. If initially set to ! NULL, it will be internally allocated and initialized to a new buffer with ! size length. Can be deallocated with FREE when no longer needed. ! ! r: on input, REAL ARRAY containing the residual vector of the iteration to ! be accelerated; on output, substituted in place with the improved residual ! vector provided by the BoostConv algorithm. ! ! length: only significant when ws points to a NULL BoostConvWS. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF ws=EMPTY BoostConvWS THEN ! allocate new workspace and exit ws=NEW BoostConvWS(length,r.LO,r.HI) WITH ws ! zero history, and set first vector equal to the current residual in(∗)=0 out(∗)=0 dd(∗)=1 rot=0 rot2=0 in(0)=r out(0)=r ELSE ! This is the actual iteration WITH ws REAL c(0..N-1) ! :projection coefficients of the current residual in(rot)=~-r ! Update in(rot), provisionally containing the residual from ! the previous iteration, to the difference between the last and the present ! residual as required by the BoostConv logic. out(rot)=~-in(rot) ! Update out(rot), containing the boosted residual ! produced by BoostConv at the previous iteration, to the difference out-in DO dd(rot,j)=in(rot)|in(j); dd(j,rot)=dd(rot,j) FOR ALL j ! Update one row and one column only of dd with the new scalar products, ! the other basis vectors being unchanged DO c(j)=in(j)|r FOR ALL j ! Project the current residual onto the non- c=PLU(dd)\~ ! orthogonal "in" basis by solving an NxN linear system IF rot=0 THEN #ifdef DiscardOldest rot=N-1 #else IF rot2=0 THEN rot2=N rot2=rot2-1 rot=rot2 #endif ELSE rot=rot-1 ! Select the next basis vector to be worked upon: circulate by decreasing ! rot by 1, and when the index hits zero reinitialize to rot2. rot2 ! in turn circulates at a lower rate. See description of the strategy. LOOP FOR ALL k in(rot,k)=r(k) r(k)=~+c∗out(∗,k) ! The "out" array actually contains out-in at this point out(rot,k)=r(k) REPEAT ! Update in place the current residual r with r+c∗(out-in), while making ! copies of the original residual into in(rot) and of the boosted residual ! into out(rot), in preparation for the next iteration. END IF END BoostConv