Commit 7877e421 authored by Cyrille Piacibello's avatar Cyrille Piacibello

So far, but not so good

parent 9c1573cb
......@@ -35,7 +35,8 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
//W and P will be stored in the same block
BlockWP<Scalar,Primary> WP(SizeBlock,B.getLeadingDim());
//Store Solution at each iteration
Block<Scalar,Primary> Y(SizeBlock,dim);
//If no Restart, Restart corresponds to the Dimension/Size of Block
HessExtended<Scalar> L(SizeBlock,B.getLeadingDim()/SizeBlock);
......@@ -63,9 +64,9 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
sizeOfEachBlock.push_back(SizeBlock);
std::vector<int> sumOfSizeOfEachBlock;
sumOfSizeOfEachBlock.push_back(SizeBlock);
int j;
//Main Loop
for(int j=1 ; j<max_iter ; ++j){
for(j=1 ; j<=max_iter ; ++j){
//Compute block Wj inside WP
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
WP.getW());
......@@ -73,7 +74,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
L.displayHessExtended("Hess Before ortho");
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim());
L.displayHessExtended("Hess After ortho");
std::cout<<"Ortho against V done"<<std::endl;
//Orthogonalisation of Wj against Pj
......@@ -81,10 +82,9 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
std::cout<<"Ortho against of W against Pj"<<std::endl;
//QR of Wj --> {Wj~,Dj}
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
std::cout<<"QR of W done"<<std::endl;
//Create Block for Y and Residual
Block<Scalar,Primary> Y(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Create Block and Residual
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
L.ComputeResidual(R1,Y,Res);
......@@ -92,6 +92,11 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//Do we really need that ?
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm();
if(MinMax.second < epsilon){
std::cout<<"Convergence ?"<<std::endl;
break;
}
//Block to Store U
Block<Scalar,Primary> U;
Res.DecompositionSVD(U,epsilon_R);
......@@ -103,7 +108,14 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
Block<Scalar,Primary> tempRpart(U.getSizeBlock(),U.getSizeBlock());
U.displayBlock("U: Kept directions");
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
if(U.getSizeBlock() == SizeBlock){//Thus no inexact breakdown occured
for(int i=0; i<SizeBlock ; ++i) {
Directions.getPtr(i)[i] = Scalar(1);
}
}else{
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
}
Directions.displayBlock("Q part of U after QR");
//Update
......@@ -112,19 +124,26 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
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} |
L.UpdateBottomLine(Directions,U.getSizeBlock());
L.displayHessExtended("hess after update bottom line");
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
//exit(0);
}
return max_iter;
//sol = Base*Y + X0
base.ComputeProduct(Y,Sol);
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
Sol.getPtr(i)[j] += X0.getPtr(i)[j];
}
}
return j;
}
#endif //BLOCKARNOLDIIB_HPP
......@@ -63,6 +63,7 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
base.addBlockDatas(V1);
......@@ -166,6 +167,7 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
globalStep++;
Hess.displayHess("Hess at the end of the loop");
}
k += nb_iter_restart;
nbRestart += 1;
......
......@@ -111,13 +111,18 @@ public:
memset(r,0,R1.getSizeBlock()*SizeBlock*(currentSize+1)*sizeof(Scalar));
int ldr = (currentSize+1)*SizeBlock;
std::cout<<"R : is "<<R1.getSizeBlock()<<" x "<<ldr<<"...\n";
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0; j<=i ; ++j){
r[ldr*i+j] = R1.getPtr(i)[j];
}
}
for(int i=0 ; i<ldr*SizeBlock ; ++i){
std::cout<<r[i]<<"\n";
}
Scalar*workingCopy = new Scalar[SizeBlock*SizeBlock*
(currentSize+1)*currentSize];
for(int i=0 ; i<currentSize ; ++i){
for(int j=0 ; j<SizeBlock ; ++j){
memcpy(&workingCopy[i*SizeBlock*ldr + j*ldr], //Dest
......@@ -126,18 +131,20 @@ public:
}
}
displayHess("Hess before gels");
callLapack_gels<Scalar>(ldr,currentSize*SizeBlock,SizeBlock,
workingCopy,ldr,r,ldr);
workingCopy,ldr,
r,ldr);
//Build return Matrix
//copy inside Yj
for(int i=0 ; i<Yj.getSizeBlock() ; ++i){
for(int j=0 ; j < currentSize*SizeBlock ; ++j){
for(int j=0 ; j < ldr ; ++j){
Yj.getPtr(i)[j] = r[ldr*i + j];
}
}
//Yj.displayBlock("Res of Gels");
Yj.displayBlock("Y resulting of Gels");
delete [] workingCopy;
delete [] r;
}
......@@ -149,11 +156,13 @@ public:
std::vector<Primary> normB){
//The Idea is to compute a matrix vector product...
Block<Scalar,Primary> blockRes(SizeBlock,(currentSize+1)*SizeBlock);
call_gemm<Scalar>((currentSize+1)*SizeBlock,SizeBlock,SizeBlock*currentSize,
datas, ldh,
Y.getPtr(), Y.getLeadingDim(),
blockRes.getPtr(),blockRes.getLeadingDim(),
1,0);
std::cout<<"About to compute Block Res:\t"<<(currentSize+1)*SizeBlock
<<"-"<<SizeBlock<<"-"<<SizeBlock*currentSize<<"\n";
call_gemm<Scalar>((currentSize+1)*SizeBlock,SizeBlock,SizeBlock*currentSize,
datas, ldh,
Y.getPtr(), Y.getLeadingDim(),
blockRes.getPtr(),blockRes.getLeadingDim(),
Scalar(1),Scalar(0));
//blockRes - R
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
......@@ -161,7 +170,7 @@ public:
blockRes.getPtr(idxRHS)[i] -= R.getPtr(idxRHS)[i];
}
}
blockRes.displayBlock("Block residual");
//Loop to get the min/max of residuals
std::pair<Primary,Primary> MinMax = std::make_pair<>(1000,0);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
......
......@@ -142,7 +142,9 @@ public:
memcpy(&r[i*ldR],gamma.getPtr(i),
sizeof(Scalar)*(i+1));
}
for(int i=0 ; i<ldR*p ; ++i){
std::cout<<r[i]<<"\n";
}
//get a working copy of F
Scalar*workingCopy = new Scalar[ldR*sumNbVectInBlck.back()];
memset(workingCopy,0,ldR*sumNbVectInBlck.back()*sizeof(Scalar));
......@@ -150,6 +152,7 @@ public:
memcpy(&workingCopy[i*ldR],&data[i*ldHess],
sizeof(Scalar)*ldR);
}
//Call least square
int gels = callLapack_gels<>(ldR,sumNbVectInBlck.back(),p,
workingCopy,ldR,
......@@ -158,17 +161,21 @@ public:
std::cout<<"Problem on Least Square\n\t gels="<<gels<<"\n";
}
//Copy back res inside Y
memcpy(outputY.getPtr(),r,sizeof(Scalar)*ldR*p);
//Call gemm : outputRes = gama - F*Y
call_gemm<>(nbLineUsed + p,p,outputY.getLeadingDim(),
for(int i=0 ; i<p ; ++i){
memcpy(outputY.getPtr(i),&r[i*ldR],sizeof(Scalar)*ldR);
}
outputY.displayBlock("Y resulting from least square");
//Call gemm : outputRes = F*Y-gama
call_gemm<>(nbLineUsed + p,p,nbVectInBlck.back(),
data,ldHess,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(-1),Scalar(0));
Scalar(1),Scalar(0));
//Add gamma to residual
for(int i=0 ; i<outputResidual.getSizeBlock();++i){
for(int j=0 ; j<i+1 ; ++j){
outputResidual.getPtr(i)[j] += gamma.getPtr(i)[j];
outputResidual.getPtr(i)[j] -= gamma.getPtr(i)[j];
}
}
outputResidual.displayBlock("Residual obtained from LeastSquare");
......
#include <random>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
#include "../src/BlockArnoldiIb.hpp"
/**
* @brief display matrix in Matlab format
* i.e [line1 ; line2 ; line3 ; ... ; line_k]
* @param : datas,
* @param m nb line
* @param n nb col
*
*/
template<class Scalar>
void display_matrix_matlab(Scalar * data,int m,int n){
std::cout<<"[ ";
for(int i=0 ; i<m ; ++i){
for(int j=0 ; j<n ; ++j){
std::cout<<" "<<data[j*m+i].real()<<"+"<<data[j*m+i].imag()<<"i ";
}
std::cout<<" ; ";
}
std::cout<<"]\n";
}
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
* useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
* size() member function in order to works with BlockGMRes algorithm
*
* An abstract class could be used to enforce that.
*/
template<class Scalar,class Primary=Scalar>
class InputMatrix{
int dim;
Scalar * data;
public:
InputMatrix(int inSize) : dim(inSize),data(nullptr){
data = new Scalar[dim*dim];
}
~InputMatrix(){
delete [] data;
}
int size(){
return dim;
}
Scalar * getPtr(){
return data;
}
void initMatrix(std::random_device& rd){
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
for(int i=0 ; i<dim ; ++i){
for(int j=0 ; j<dim ; ++j){
data[dim*i+j] = {dis(gen),dis(gen)};
}
}
}
void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out){
call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
data,dim,
in.getPtr(), in.getLeadingDim(),
out.getPtr(),out.getLeadingDim(),
Scalar(1),Scalar(0));
}
void MatBaseProduct(Scalar * toRead, int nbRHS, Scalar* ptrToWrite){
call_gemm<Scalar>(dim,nbRHS, dim,
data,dim,
toRead,dim,
ptrToWrite,dim,
Scalar(1),Scalar(0));
}
bool useRightPreCond(){
return false;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
for(int i=0 ; i<dim ; ++i){
output.getPtr(j)[i] = input.getPtr(j)[i]*(Scalar(1)/data[i*dim+i]);
}
}
}
void preCondBaseProduct(Scalar * ptr,Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
for(int i=0 ; i<dim ; ++i){
output.getPtr(j)[i] = ptr[j*dim+i]*(Scalar(1)/(data[i*dim+i]));
}
}
}
};
/**
* @brief This executable uses BlockArnoldiIb procedure to generate a
* base and compute a solution to a system.
*
*/
int main(int ac, char ** av){
std::vector<std::string > args(av,av+ac);
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock=4,nbBlock=5;
if (args.size()>1){
SizeBlock = std::stoi(args[1]);
}
if (args.size()>2){
nbBlock = std::stoi(args[2]);
}
int dim = nbBlock*SizeBlock;
using Primary=double;
using Scalar=std::complex<double>;
//Random engine used to generate matrix
std::random_device rd;
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
//Create matrix
InputMatrix<Scalar,Primary> mat(dim);
mat.initMatrix(rd);
//Create X_Exact
Block<Scalar,Primary> X_Exact(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
X_Exact.getPtr(i)[j] = {dis(gen),dis(gen)};
}
}
//Compute RHS as Mat*X_Exact
Block<Scalar,Primary> RHS(SizeBlock,dim);
mat.MatBlockVect(X_Exact,RHS);
//Compute A starting Block
Block<Scalar,Primary> X0(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock-1 ; ++i){ //last col not filled
for(int j=0 ;j<dim; ++j){
X0.getPtr(i)[j] = {dis(gen),dis(gen)};
}
}
//fill last col with one vect from solution
int last_indice = SizeBlock-1;
memcpy(X0.getPtr(last_indice),X_Exact.getPtr(0),
sizeof(Scalar)*X0.getLeadingDim());
Block<Scalar,Primary> sol(SizeBlock,dim);
Block<Scalar,Primary> sol_ref(SizeBlock,dim);
display_matrix_matlab(mat.getPtr(),dim,dim);
std::cout<<"\n\n";
display_matrix_matlab(RHS.getPtr(),dim,SizeBlock);
std::cout<<"\n\n";
display_matrix_matlab(X0.getPtr(),dim,SizeBlock);
std::cout<<"\n\n";
Logger<Primary> logger;
int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
X0,5,sol);
std::cout<<"Return value : "<<nb<<"\n";
//Check solution
{
//Compute what i need for my metrics
//MatxSol : Mat*Sol_computed
Block<Scalar,Primary> MatxSol(SizeBlock,dim);
mat.MatBlockVect(sol,MatxSol);
//MatxSolMinusB : Mat*Sol_computed-B
Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
//X_Diff : X_Exact - Sol_computed
Block<Scalar,Primary> X_Diff(SizeBlock,dim);
MatxSolMinusB.CopyBlock(MatxSol);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
for(int j=0 ; j<dim; ++j){
MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
sol.getPtr(idxRHS)[j];
}
}
//Then display for each vect
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
std::cout<<"Metrics on "<<idxRHS<<"\n";
Primary ForwardError = X_Diff.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error (on X) \t" << ForwardError <<"\n";
Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Second Member Error \t" << SecondMemberError <<"\n";
Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
Primary normB = RHS.getNorm(idxRHS);
Primary NormwiseBackwardError = SecondMemberError / normB;
std::cout <<"["<<idxRHS<<"]\t"
<<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
}
}
int nb2 = BlockGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,5,
sol_ref,5,logger);
std::cout<<"Return value : "<<nb<<"\n";
//Check solution
{
//Compute what i need for my metrics
//MatxSol : Mat*Sol_computed
Block<Scalar,Primary> MatxSol(SizeBlock,dim);
mat.MatBlockVect(sol_ref,MatxSol);
//MatxSolMinusB : Mat*Sol_computed-B
Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
//X_Diff : X_Exact - Sol_computed
Block<Scalar,Primary> X_Diff(SizeBlock,dim);
MatxSolMinusB.CopyBlock(MatxSol);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
for(int j=0 ; j<dim; ++j){
MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
sol_ref.getPtr(idxRHS)[j];
}
}
//Then display for each vect
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
std::cout<<"Metrics on "<<idxRHS<<"\n";
Primary ForwardError = X_Diff.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error (on X) \t" << ForwardError <<"\n";
Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Second Member Error \t" << SecondMemberError <<"\n";
Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
Primary normB = RHS.getNorm(idxRHS);
Primary NormwiseBackwardError = SecondMemberError / normB;
std::cout <<"["<<idxRHS<<"]\t"
<<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
}
}
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment