BlockArnoldiIb.hpp 6.55 KB
Newer Older
1 2 3
#ifndef BLOCKARNOLDIIB_HPP
#define BLOCKARNOLDIIB_HPP

4 5
#include "BlockWP.hpp"
#include "HessExtended.hpp"
6 7
#include "Logger.hpp"
#include "Base.hpp"
8

9 10 11 12 13
/**
 * @brief : Arnoldi iteration with inexact breakdowns
 *
 * Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
 */
14
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
15
int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
16
                   int max_iter, Block<Scalar,Primary>& Sol,Logger<Primary,N>& log,
17
                   Primary epsilon = 0.0000001){
18 19
    double tic; //main tic
    double local_tic;
20 21
    //For monitoring purpose
    int nbDeflDir = 0;
22

Cyrille Piacibello's avatar
Cyrille Piacibello committed
23
    int SizeBlock = B.getSizeBlock();
24
    int dim = A.size();
25

26 27
    //store norms
    std::vector<Primary> normB;
28
    for(int i=0 ; i< SizeBlock ; ++i){
29
        normB.push_back(B.getNorm(i));
30 31 32
    }
    //Compute epsilon_R (step 1)
    Primary epsilon_R = normB.front();
33
    for(auto &a : normB){
34 35 36 37 38
        if(a<epsilon_R){
            epsilon_R = a;
        }
    }
    epsilon_R *= epsilon;
Cyrille Piacibello's avatar
Cyrille Piacibello committed
39 40

    //Here will the restart loop
41
    Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
42 43
    //W and P will be stored in the same block
    BlockWP<Scalar,Primary> WP(SizeBlock,B.getLeadingDim());
Cyrille Piacibello's avatar
Cyrille Piacibello committed
44
    //Store Solution at each iteration
45
    Block<Scalar,Primary> Y(SizeBlock,dim+SizeBlock);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
46 47
    //If no Restart, Restart corresponds to the Dimension/Size of Block
    HessExtended<Scalar> L(SizeBlock,B.getLeadingDim()/SizeBlock);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
48

49 50
    log.init_log(SizeBlock);

51
    //Compute first block residual (step 1)
52
    Block<Scalar,Primary> R0(SizeBlock,dim);
53 54 55
    //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
56 57
        for(int k=0 ; k<dim ; ++k){
            R0.getPtr(i)[k] = B.getPtr(i)[k] - R0.getPtr(i)[k];
58 59 60 61
        }
    }
    //V1 will be first block of base, and R1, will be used for least
    //square pb, in criteria verification
62 63
    Block<Scalar,Primary> V1(SizeBlock,dim);
    Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
64 65
    //Q-R Decomposition of (R0)
    R0.ComputeQRFacto(V1,R1);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
66 67
    //Fisrt block of base added
    base.addBlockDatas(V1);
68

Cyrille Piacibello's avatar
Cyrille Piacibello committed
69
    //We store the size of each block -- Could be delegated to a class ?
Cyrille Piacibello's avatar
Cyrille Piacibello committed
70 71
    std::vector<int> sizeOfEachBlock;
    sizeOfEachBlock.reserve(max_iter);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
72 73 74
    sizeOfEachBlock.push_back(SizeBlock);
    std::vector<int> sumOfSizeOfEachBlock;
    sumOfSizeOfEachBlock.push_back(SizeBlock);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
75
    int j;
76 77 78 79
    //Add first ite
    std::pair<Primary,Primary> MinMax = R0.getMinMaxNorm(normB);
    log.add_ite(-1,MinMax.first,MinMax.second,B.getSizeBlock());

Cyrille Piacibello's avatar
Cyrille Piacibello committed
80
    //Main Loop
81 82
    for(j=0 ; j<max_iter ; ++j){
        tic=Logger<Primary,N>::GetTime();
83 84
        std::cout<<"################## Start of Iteration #################\n";

85
        //Compute block Wj inside WP
86
        local_tic=Logger<Primary,N>::GetTime();
87 88
        A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
                         WP.getW());
89
        log.add_step(1,local_tic);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
90
        //Ortho against V and store coeff inside new block column of Hess L
91
        //L.displayHessExtended("Hess Before ortho");
92
        Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
93
        local_tic=Logger<Primary,N>::GetTime();
94
        base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
95
        log.add_step(2,local_tic);
96
        std::cout<<"Ortho against V and P"<<std::endl;
Cyrille Piacibello's avatar
Cyrille Piacibello committed
97

Cyrille Piacibello's avatar
Cyrille Piacibello committed
98 99
        //Orthogonalisation of Wj against Pj
        WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
100

Cyrille Piacibello's avatar
Cyrille Piacibello committed
101 102
        //QR of Wj --> {Wj~,Dj}
        WP.QRofW(L.getPtrToD(),L.getLeadingDim());
103
        //L.template displayHessExtendedBitMap<Primary>("L after Orthogonalization");
104
        std::cout<<"QR of W done"<<std::endl;
Cyrille Piacibello's avatar
Cyrille Piacibello committed
105 106

        //Create Block and Residual
Cyrille Piacibello's avatar
Cyrille Piacibello committed
107 108
        Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
        //Compute residual
109
        local_tic=Logger<Primary,N>::GetTime();
Cyrille Piacibello's avatar
Cyrille Piacibello committed
110
        L.ComputeResidual(R1,Y,Res);
111
        log.add_step(3,local_tic);
112 113
        std::cout<<"\nResidual Computed ...\n";

Cyrille Piacibello's avatar
Cyrille Piacibello committed
114
        //Do we really need that ?
115
        MinMax = Res.getMinMaxNorm(normB);
116 117 118
        std::cout<<"Ite ["<<j<<"] Done\t"<<
            "Min\t"<<MinMax.first<<"\t"<<
            "Max\t"<<MinMax.second<<"\n";
119
        log.add_ite(j,MinMax.first,MinMax.second,WP.getSizeW());
Cyrille Piacibello's avatar
Cyrille Piacibello committed
120

Cyrille Piacibello's avatar
Cyrille Piacibello committed
121
        if(MinMax.second < epsilon){
122
            std::cout<<"Convergence !!"<<std::endl;
Cyrille Piacibello's avatar
Cyrille Piacibello committed
123 124 125
            break;
        }

126

Cyrille Piacibello's avatar
Cyrille Piacibello committed
127 128
        //Block to Store U
        Block<Scalar,Primary> U;
129
        local_tic=Logger<Primary,N>::GetTime();
Cyrille Piacibello's avatar
Cyrille Piacibello committed
130
        Res.DecompositionSVD(U,epsilon_R);
131
        log.add_step(4,local_tic);
132 133
        std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
                 <<U.getSizeBlock()<<" !"<<std::endl;
Cyrille Piacibello's avatar
Cyrille Piacibello committed
134

135

136 137
        //Here : If no inexact breakdowns occured, updates procedures are easiers
        if(U.getSizeBlock() == SizeBlock){
138
            //Augment base : [P_{j-1},~Wj] = [~Wj]
139
            base.addBlockDatas(WP);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
140
        }else{
141 142 143 144
            if(U.getSizeBlock() + nbDeflDir != SizeBlock){
                nbDeflDir++;
            }
            std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
145 146
            //Create Block to store \W
            Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
147
            Block<Scalar,Primary> tempRpart;
148 149

            //Qr facto of bottom block of U
Cyrille Piacibello's avatar
Cyrille Piacibello committed
150
            U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
151 152 153 154 155 156 157 158 159 160

            //Update
            //V_{j+1} = [P_{j-1},~Wj] * Directions_1
            WP.ComputeProduct(Directions,U.getSizeBlock(),
                              base.getLastColPtr(U.getSizeBlock()),base.getLeadingDim());
            //Pj = [P_{j-1},~Wj] * Directions_2
            WP.ComputeProduct(Directions,U.getSizeBlock());
            //L.displayHessExtended("hess before update bottom line");
            //L_{j+1,} = | Directions_1^{H} | * H_j
            //Gj =       | Directions_2^{H} |
161 162
            L.UpdateBottomLine(Directions);
            //L.template displayHessExtendedBitMap<Primary>("L after update");
Cyrille Piacibello's avatar
Cyrille Piacibello committed
163
        }
Cyrille Piacibello's avatar
Cyrille Piacibello committed
164 165
        sizeOfEachBlock.push_back(U.getSizeBlock());
        sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
166

167 168
        std::cout<<"################## End of Iteration ##################\n"
                 <<"Base contains exactly :: "<<base.getNbVectUsed()<<"\t!!\n";
169
        log.add_step(0,tic);
Cyrille Piacibello's avatar
Cyrille Piacibello committed
170
    }
Cyrille Piacibello's avatar
Cyrille Piacibello committed
171 172 173
    //sol = Base*Y + X0
    base.ComputeProduct(Y,Sol);
    for(int i=0 ; i<SizeBlock ; ++i){
174 175
        for(int k=0 ; k<dim ; ++k){
            Sol.getPtr(i)[k] += X0.getPtr(i)[k];
Cyrille Piacibello's avatar
Cyrille Piacibello committed
176 177 178 179
        }
    }

    return j;
180 181 182
}

#endif //BLOCKARNOLDIIB_HPP