BlockGMRes.hpp 8.67 KB
Newer Older
1 2 3 4 5
#ifndef BLOCKGMRES_HPP
#define BLOCKGMRES_HPP


#include "LapackInterface.hpp"
6
#include "HessStandard.hpp"
7 8
#include "Base.hpp"
#include "Block.hpp"
9
#include "Logger.hpp"
10 11 12 13 14 15 16

/**
 * @brief Computation of the GMRES Block.
 *
 * Ref : Saad P.241
 * Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
 */
17
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
18
int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
19
               //Matrix preCond Deal with that later
20
               int max_iter, Block<Scalar,Primary>& Sol, int nb_iter_restart,
21
               Logger<Primary,5>& log,
22 23
               Primary epsilon = 0.00001){

24 25 26 27 28
    if(A.useRightPreCond()){
        std::cout<<"Right Pre Conditionner used\n";
    }else{
        std::cout<<"NO - Right Pre Conditionner used\n";
    }
29 30 31 32 33 34 35 36

    //Global variables
    int dim = A.size();
    int SizeBlock = B.getSizeBlock();

    //store norms
    std::vector<Primary> normB;
    for(int i=0 ; i< SizeBlock ; ++i){
37
        normB.push_back(B.getNorm(i));
38 39 40
    }

    //Residual vectors needed for computing the solution
41
    Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
42 43
    //Add a loop for restart
    int k=0,nbRestart=0, globalStep=0;
44
    log.init_log(SizeBlock);
45 46 47
    while(k<max_iter){
        std::cout<<"\t\t\tRestarting number "<<nbRestart<<" ("<<nb_iter_restart<<")\n";

48 49
        //        X0.displayBlock("New X0");

50
        Block<Scalar,Primary> R0(SizeBlock,dim);
51 52 53 54 55 56 57 58 59
        //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 j=0 ; j<dim ; ++j){
                R0.getPtr(i)[j] = B.getPtr(i)[j] - R0.getPtr(i)[j];
            }
        }

        //Storage needed for computating
60
        Base<Scalar,Primary> base(SizeBlock,nb_iter_restart+1,dim);
61 62 63
        Hessenberg<Scalar> Hess(nb_iter_restart,SizeBlock);
        //V1 will be first block of base, and R1, will be used for least
        //square pb, in criteria verification
64 65
        Block<Scalar,Primary> V1(SizeBlock,dim);
        Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
66 67
        //Q-R Decomposition of (R0)
        R0.ComputeQRFacto(V1,R1);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
68

69 70 71 72 73 74
        //Fisrt block of base added
        base.addBlockDatas(V1);

        std::cout<<"Starting main loop"<<std::endl;
        for(int j=0 ; (j<nb_iter_restart) && (j+k)<= max_iter ; ++j){
            //Compute block Wj
75
            Block<Scalar,Primary> Wj(SizeBlock,dim);
76

77
            if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
78
                Block<Scalar,Primary> temp(SizeBlock,dim);
79 80 81
                A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
                A.MatBlockVect(temp,Wj);
            }else{//Wj = A*V[j] without preCond
82 83
                A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),
                                 Wj.getSizeBlock(),Wj.getPtr());
84
            }
85 86 87 88 89 90 91 92

            //Where to write in hess
            Scalar * toWrite = Hess.getPtr(j);
            base.ComputeMGSOrtho(Wj,toWrite,Hess.getLeadingDim());

            std::cout<<"Ortho done"<<std::endl;

            //QR facto of orthogonalized Wj
93 94
            Block<Scalar,Primary> Q(SizeBlock,dim);
            Block<Scalar,Primary> R(SizeBlock,SizeBlock);
95 96 97
            Wj.ComputeQRFacto(Q,R);
            //Copy Block R inside Hess
            Hess.addRPartToHess(R);
98

99 100 101 102
            std::cout<<"QR done"<<std::endl;

            //Tests over criterion
            Hess.ComputeBlockResidual(R1,Yj);
103 104

            //Here we compute Base*Y +X0 and display res
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
            // {
            //     //if no precond
            //     //    X0 = X0 + Base*Yj
            //     Block<Scalar,Primary> Resultat{SizeBlock,dim};
            //     //Yj.displayBlock("Yj at current iter");
            //     base.ComputeProduct(Yj,Resultat);
            //     for(int i=0 ; i<SizeBlock ; ++i){
            //         for(int j=0 ; j<dim ; ++j){
            //             Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
            //         }
            //     }
            //     std::cout<<"#########\t"<<j<<"\t#########\n";
            //     //Resultat.displayBlock("Current Xm");

            //     //Then we compute real residual
            //     Block<Scalar,Primary> AxSol(SizeBlock,dim);
            //     A.MatBlockVect(Resultat,AxSol);
            //     for(int i=0 ; i<SizeBlock ; ++i){
            //         for(int j=0 ; j<dim ; ++j){
            //             AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
            //         }
            //     }
            //     for(int i=0 ; i<SizeBlock ; ++i){
            //         std::cout<<"PREFIX "<<i<<"\t"<<j<<"\t"<<AxSol.getNorm(i)<<"\n";
            //     }
            // }
131 132


133 134 135 136
            std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;

            {//Display residual to see convergence
                std::pair<Primary,Primary> MinMax
137
                    = Hess.template displayBlockResidual<Primary>(Yj,R1,normB);
138
                log.add_ite(j,MinMax.first,MinMax.second);
139 140 141 142 143 144
                std::cout<<"["<<j<<"]\t"<<"Min "<<MinMax.first
                         <<"\tMax "<<MinMax.second<<"\n";
                if(MinMax.second <= epsilon){
                    //Compare to the real residual
                    std::pair<Primary,Primary> RealMinMax
                        = std::make_pair<Primary,Primary>(100,0);
145 146 147

                    //Generate block Sol = Base*Yj
                    //if Precond : block Sol = preCond*Base*yj
148
                    Block<Scalar,Primary> blockSol(SizeBlock,dim);
149
                    if(A.useRightPreCond()){
150
                        Block<Scalar,Primary> temp(SizeBlock,dim);
151 152 153 154 155
                        base.ComputeProduct(Yj,temp);
                        A.preCondBlockVect(temp,blockSol);
                    }else{
                        base.ComputeProduct(Yj,blockSol);
                    }
156 157

                    //Compute addition with X0
158 159 160
                    for(int i=0 ; i<SizeBlock ; ++i){
                        for(int j=0 ; j<dim ; ++j){
                            blockSol.getPtr(i)[j] += X0.getPtr(i)[j];
161 162 163
                        }
                    }
                    //Compute A*sol computed
164
                    Block<Scalar,Primary> AxSol(SizeBlock,dim);
165 166 167 168 169
                    A.MatBlockVect(blockSol,AxSol);

                    //Compute real residual
                    for(int i=0 ; i<SizeBlock ; ++i){
                        for(int j=0 ; j<dim ; ++j){
170
                            AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
171 172 173
                        }
                    }
                    //For each vector, compute norm
174
                    for(int i = 0 ; i<SizeBlock ; ++i){
175
                        auto norm = AxSol.getNorm(i)/normB[i];
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
                        if(norm > RealMinMax.second){
                            RealMinMax.second = norm;
                        }
                        if(norm < RealMinMax.first){
                            RealMinMax.first = norm;
                        }
                    }//If convergence
                    if(RealMinMax.second <= epsilon){
                        std::cout<<"Convergence !\t Step : "<<j+k<<std::endl;
                        std::cout<<"Min "<<RealMinMax.first
                                 <<"\tMax "<<RealMinMax.second<<"\n";
                        //Modify k in order to stop all iterations
                        k = max_iter;
                        //Base augmentation
                        base.addBlockDatas(Q);
                        break;
                    }else{
                        std::cout
                            <<"Residual is small enough, but not the real's one \nMin2"
                            <<RealMinMax.first<<"\tMax2"<<RealMinMax.second<<"\n";
                    }
                }
            }
            //Base augmentation
            base.addBlockDatas(Q);
            std::cout<<"Base has been augmented\n";
            globalStep++;
        }
        k += nb_iter_restart;
        nbRestart += 1;
        //Compute new starting block X0
207 208 209 210
        //if no precond
        //    X0 = X0 + Base*Yj
        //else
        //    X0 = X0 + preCond*Base*yj
211
        Block<Scalar,Primary> newX0(SizeBlock,dim);
212
        if(A.useRightPreCond()){
213
            Block<Scalar,Primary> temp(SizeBlock,dim);
214 215 216 217 218
            base.ComputeProduct(Yj,temp);
            A.preCondBlockVect(temp,newX0);
        }else{
            base.ComputeProduct(Yj,newX0);
        }
219
        //Copy back
220 221
        //End of use of Yj, reset it
        Yj.reset();
222 223 224 225 226 227
        for(int i=0 ; i<SizeBlock ; ++i){
            for(int j=0 ; j<dim ; ++j){
                X0.getPtr(i)[j] = X0.getPtr(i)[j] + newX0.getPtr(i)[j];
            }
        }
    }
228
    Sol.CopyBlock(X0);
229 230 231
    return globalStep;
}

232
#endif //BLOCKGMRES_HPP