Arnoldi_IB.hpp 10.6 KB
Newer Older
1 2 3 4 5 6 7
#ifndef ARNOLDI_IB_HPP
#define ARNOLDI_IB_HPP

#include "Base.hpp"
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Utils.hpp"
8
#include "Logger.hpp"
9 10 11 12 13 14 15 16 17

/**
 * @brief : Arnoldi iterations with inexact breakdowns
 *
 * @param A : User Matrix
 * @param X0 : Needed for computing real solution if convergence
 * @param B : Needed for computing real residual if convergence
 * @param max_ite : number of iteration to do (at max)
 * @param epsilon : tolerance for Residual
18 19
 * @param log : logger
 *
20 21 22 23 24 25 26 27
 * @return Number of iterations done
 *
 * Solution will be stored inside X0 at the end of computation.
 *
 * Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
 */
template<class Matrix,class Scalar,class Primary=Scalar>
struct Arnoldi_IB{
28 29 30 31
    static ArnReturn Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
                                    Block<Scalar,Primary>& B,
                                    int max_iter,int MaxKSize, Primary epsilon,
                                    Logger<Primary>& log){
32 33 34 35 36
        std::cout<<"########### Arnoldi with Inexact Breakdown ############\n";
        //Initial Size of block
        int SizeBlock = B.getSizeBlock();
        //Size of matrix
        int dim = A.size();
37 38 39

        //Return struct
        ArnReturn returned_struct;
40
        returned_struct.sizeKSpace = 0;
41 42

        //Set the hessenberg
43
        HessExtended<Scalar,Primary> L(MaxKSize,SizeBlock);
44 45 46 47 48 49
        Base<Scalar,Primary> base(SizeBlock,MaxKSize,dim);
        //Create Block to Store W and P (if IB happens)
        BlockWP<Scalar,Primary> WP(SizeBlock,dim);
        //Create Block to Store Y
        Block<Scalar,Primary> Y(SizeBlock,dim+SizeBlock);

50 51 52 53 54 55 56 57 58
        //Compute first block residual (step 1)
        Block<Scalar,Primary> R0(SizeBlock,dim);
        //Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
        A.MatBlockVect(X0,R0);             //R0 = A*X0
        for(int i=0 ; i<SizeBlock ; ++i){  //R0 = B - R0
            for(int k=0 ; k<dim ; ++k){
                R0.getPtr(i)[k] = B.getPtr(i)[k] - R0.getPtr(i)[k];
            }
        }
59 60
        Block<Scalar,Primary> q0(SizeBlock,dim);
        Block<Scalar,Primary> r0(SizeBlock,SizeBlock);
61
        //Q-R Decomposition of (R0)
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
        R0.ComputeQRFacto(q0,r0);

        //here, insert inexact breakdown on R0
        {
            //Copy r0 because needed to compute Lambda, but is
            //destroyed if SVD is called
            Block<Scalar,Primary> r0_cpy(SizeBlock,SizeBlock);
            r0_cpy.CopyBlock(r0);

            //U is of full size, in order to get all the left singular vectors
            Block<Scalar,Primary> U{r0.getSizeBlock(),r0.getLeadingDim()};
            //SVD
            int nbKDir = r0_cpy.DecompositionSVD(U,epsilon);
            //Test if needed to launch algorithm again
            if(nbKDir == 0){
                std::cout<<"\t\tINEXACT BREAKDOWN on R0\n";
                std::cout<<"\t\tNb Kept Dir :"<<nbKDir<<"/"<<SizeBlock<<std::endl;
                returned_struct.sizeKSpace = 0;
                returned_struct.nbIteDone  = 0;
                returned_struct.hasConverged  = true;
                return returned_struct;
            }

            if(nbKDir == SizeBlock){
                std::cout<<"\t\tNo INEXACT BREAKDOWN on R0\n";
                //First block of base added
                base.addBlockDatas(q0);
                L.initLambda(r0);
            }else{
                std::cout<<"\t\tINEXACT BREAKDOWN on R0\n";
                std::cout<<"\t\tNb Kept Dir :"<<nbKDir<<"/"<<SizeBlock<<std::endl;
                //Block empty that will contains the first base vectors
                Block<Scalar,Primary> V1ToBeAdded{nbKDir,dim};

                //V1ToBeAdded = V1 * U1
                call_gemm<>(dim,nbKDir,SizeBlock,
                            q0.getPtr(),q0.getLeadingDim(),
                            U.getPtr(),U.getLeadingDim(),
                            V1ToBeAdded.getPtr(),V1ToBeAdded.getLeadingDim(),
                            Scalar(1),Scalar(0));

                //P = V1*U2
                call_gemm<>(dim,SizeBlock-nbKDir,SizeBlock,
                            q0.getPtr(),q0.getLeadingDim(),
                            U.getPtr(nbKDir),U.getLeadingDim(),
                            WP.getP(),WP.getLeadingDim(),
                            Scalar(1),Scalar(0));
                WP.incSizeP(SizeBlock-nbKDir);

                base.addBlockDatas(V1ToBeAdded);

                //Create a Lambda1 = [U1,U2]^{H} * R0;
                Block<Scalar,Primary> Lambda1{SizeBlock,SizeBlock};

                call_Tgemm<>(SizeBlock,SizeBlock,SizeBlock,
                             U.getPtr(),U.getLeadingDim(),
                             r0.getPtr(),r0.getLeadingDim(),
                             Lambda1.getPtr(),Lambda1.getLeadingDim(),
                             Scalar(1),Scalar(0));

                L.initLambda(Lambda1);

                //Create and initialize Phi (\Phi_1)
                L.initPhi(nbKDir);
            }
        }

        //We store the size of each block -- Could be delegated to a class ?
        std::vector<int> sumOfSizeOfEachBlock;
        sumOfSizeOfEachBlock.push_back(SizeBlock - WP.getSizeP());
132 133

        //Nb directions discarded
134 135
        int nbDeflDir = WP.getSizeP();
        bool convergence = (nbDeflDir == SizeBlock);
136 137
        //Counter
        int j=0;
138

139 140
        while(j<max_iter){        //Main Loop
            double timestamp = Logger<Primary>::GetTime();
141
            std::cout<<"\t\tStart of Iteration\n";
142 143 144 145 146 147 148

            //Get new starting point, and test if restart is triggered
            Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
            if(nullptr == toWrite){ //No more room for this iteration
                break;
            }

149 150 151 152 153 154 155 156
            if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
                Block<Scalar,Primary> temp(SizeBlock-nbDeflDir,dim);
                A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
                A.MatBlockVect(temp,WP,WP.getSizeP());
            }else{//Wj = A*V[j] without preCond
                A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
                                 WP.getW());
            }
157
            returned_struct.sizeKSpace += WP.getSizeW();
158 159 160 161 162 163 164 165 166 167 168 169
            //Ortho against V and store coeff inside new block column of Hess L
            base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
            std::cout<<"Ortho against V and P"<<std::endl;

            //Orthogonalisation of Wj against Pj
            WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
            //QR of Wj --> {Wj~,Dj}
            WP.QRofW(L.getPtrToD(),L.getLeadingDim());
            std::cout<<"QR of W done"<<std::endl;

            //Create Block and Residual
            Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
170

171
            //Compute residual
172
            L.ComputeResidual(Y,Res);
173 174
            std::cout<<"\nResidual Computed ...\n";

175 176
            Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);

177 178 179 180 181 182 183 184
            //Do we really need that ?
            std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm();
            std::cout<<"Ite ["<<j<<"] Done\t"<<
                "Min\t"<<MinMax.first<<"\t"<<
                "Max\t"<<MinMax.second<<"\n";
            if(MinMax.second < epsilon){
                //here, add a test on the real residual
                if(CheckCurrentSolution(A,B,X0,base,Y,epsilon)){
185 186 187
                    //Adding last ite
                    log.add_ite(j,WP.getSizeW(),MinMax.first,MinMax.second,timestamp);
                    convergence = true;
188 189 190 191 192 193
                    break;
                }
            }
            //Block to Store U
            Block<Scalar,Primary> U;

194
            int nbDirKept = Res.DecompositionSVD(U,epsilon);
195
            std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
196
                 <<nbDirKept<<" !"<<std::endl;
197

198
            //Here : If no inexact breakdowns occured, updates procedures are easier
199 200 201 202 203 204 205 206 207
            if(U.getSizeBlock() == SizeBlock){
                //Augment base : [P_{j-1},~Wj] = [~Wj]
                base.addBlockDatas(WP);
            }else{
                if(U.getSizeBlock() + nbDeflDir != SizeBlock){
                    nbDeflDir = SizeBlock-U.getSizeBlock();
                }
                std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
                //Create Block to store \W
208

209 210 211
                Block<Scalar,Primary> tempRpart;

                //Qr facto of bottom block of U
212 213 214
                std::cout<<"About to compute QR of bottom part of U\n"
                         <<"Starting point : "<<sumOfSizeOfEachBlock.back()
                         <<"\n";
215
                U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
216 217
                //Update Phi with directions (W1W2), doc Annexe Eq 2.7
                L.updatePhi(Directions,nbDirKept);
218 219 220

                //Update
                //V_{j+1} = [P_{j-1},~Wj] * Directions_1
221 222 223 224 225 226 227 228
                Scalar * baseToWrite = base.getLastColPtr(U.getSizeBlock());
                if(baseToWrite){
                    WP.ComputeProduct(Directions,U.getSizeBlock(),
                                      baseToWrite,
                                      base.getLeadingDim());
                }else{
                    break;
                }
229 230 231 232 233 234
                //Pj = [P_{j-1},~Wj] * Directions_2
                WP.ComputeProduct(Directions,U.getSizeBlock());
                //L_{j+1,} = | Directions_1^{H} | * H_j
                //Gj =       | Directions_2^{H} |
                L.UpdateBottomLine(Directions);
            }
235
            //sizeOfEachBlock.push_back(U.getSizeBlock());
236 237 238
            sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());

            std::cout<<"################## End of Iteration ##################\n\n";
239 240
            log.add_ite(j,WP.getSizeW(),MinMax.first,MinMax.second,timestamp);
            ++j;
241
        }
242
        if(!convergence){ //No convergence, we need to write current sol into X0
243 244 245 246 247 248 249 250 251 252 253 254 255 256
            Block<Scalar,Primary> Sol{SizeBlock,dim};
            if(A.useRightPreCond()){
                Block<Scalar,Primary> temp(SizeBlock,dim);
                base.ComputeProduct(Y,temp);
                A.preCondBlockVect(temp,Sol);
            }else{//No pre Cond used
                base.ComputeProduct(Y,Sol);
            }
            for(int i=0 ; i<SizeBlock ; ++i){
                for(int k=0 ; k<dim ; ++k){
                    X0.getPtr(i)[k] += Sol.getPtr(i)[k];
                }
            }
        }
257 258 259 260
        std::cout<<"################# Iterations done ... "<<j<<"\n";
        returned_struct.nbIteDone = j;
        returned_struct.hasConverged = convergence;
        return returned_struct;
261 262 263 264 265 266
    }
};



#endif //ARNOLDI_IB_HPP