Commit b3976813 authored by Cyrille Piacibello's avatar Cyrille Piacibello

Inexact Breakdown in progress, member func in BlockWP complete, HessExtended almost finished

parent a9ef2a57
......@@ -77,8 +77,8 @@ public:
memcpy(R.getPtr(i),&datas[i*ldb],(i+1)*sizeof(Scalar));
}
//Call to orgqr //orgqr does not exist in Lapack for complex,
//but Lapack has a counter part called ungqr, which is called
//if Scalar is a complex type
//but Lapack has a complex counter part called ungqr, which is
//called if Scalar is a complex type
int orgqr = callLapack_orgqr(ldb,SizeBlock,datas,ldb,tau);
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
......@@ -107,13 +107,13 @@ public:
/**
* @brief Copy input Block datas inside
*/
void CopyBlock(Block<Scalar,Primary>& in, int nbColToCopy,int start=0){
void CopyBlock(Block<Scalar,Primary>& in){
//Check
if(in.getLeadingDim() != ldb){
std::cout<<"Leading Dim different between input block and current block\n";
exit(0);
}
if(nbColToCopy != SizeBlock){
if(in.getSizeBlock() != SizeBlock){
std::cout<<"SizeBlock different between input clo to copy"
<<" and current block\n";
exit(0);
......@@ -138,7 +138,7 @@ public:
//What I need:
char jobu = 'S'; //jobu : S --> Only first size_block vectors computed
char jobvt = 'N'; //jobvt : N --> no part of V^{H} is computed
Scalar* Vt; //Won't be used
Scalar* Vt = nullptr; //Won't be used
int ldVt = 1; //Won't be used
//Needed in case Lapack did not converge
Primary * superb = new Primary[output.getSizeBlock()-1];
......
......@@ -188,7 +188,7 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
}
}
}
Sol.CopyBlock(X0,X0.getSizeBlock());
Sol.CopyBlock(X0);
return globalStep;
}
......
......@@ -61,6 +61,26 @@ public:
* Q. R part generated will go inside input toWrite ptr
*/
void QRofW(Scalar * toWriteR, int ldR){
//Create temporary tau
Scalar * tau = new Scalar[getSizeW()];
memset(tau,0,sizeof(Scalar)*getSizeW());
//3 steps::
//1 - Call qr facto on W bloc
int geqrf = callLapack_geqrf<>(Parent::getLeadingDim(),getSizeW(),getW(),
Parent::getLeadingDim(),tau);
if(geqrf != 0){
std::cout<<"Problem on QR facto\n";
}
//2 - get the R part and write it in input
for(int i=0 ; i<getSizeW() ; ++i){
memcpy(&toWriteR[i*ldR],&getW()[i*Parent::getLeadingDim()],
(i+1)*sizeof(Scalar));
}
//3 - Call orgqr on W block
int orgqr = callLapack_orgqr(Parent::getLeadingDim(),getSizeW(),
getW(),Parent::getLeadingDim(),tau);
delete [] tau;
}
/**
......
#ifndef HESSEXTENDED_HPP
#define HESSEXTENDED_HPP
#include "Block.hpp"
#include "BlockWP.hpp"
#include "LapackInterface.hpp"
/**
......@@ -36,12 +36,12 @@ class HessExtended{
public:
HessExtended(int inP,int inRestart) : p(inP),
nbTotalCol(p*restart),
ldHess(p*(restart+1)),
nbTotalCol(p*inRestart),
ldHess(p*(inRestart+1)),
nbBlckColUsed(0),
nbLineUsed(0),data(nullptr){
data = new Scalar[nbTotalLine*nbTotalCol];
memset(data,0,sizeof(Scalar)*nbTotalLine*nbTotalCol);
data = new Scalar[ldHess*nbTotalCol];
memset(data,0,sizeof(Scalar)*ldHess*nbTotalCol);
}
~HessExtended(){
......@@ -59,7 +59,8 @@ public:
Scalar * getNewStartingCol(int inSizeBlock){
int ref = sumNbVectInBlck.back();
nbVectInBlck.push_back(inSizeBlock);
sumNbVectInBlck.push_back(ref.inSizeBlock);
sumNbVectInBlck.push_back(ref+inSizeBlock);
nbBlckColUsed++;
return &data[ref*ldHess];
}
......@@ -67,24 +68,41 @@ public:
* @brief Remove the last part (containing H), and filling it with
* W^{H}*H
*
* @param Wj
* @param W1 Block
*/
template<typename Primary>
void UpdateBottomLine(Block<Scalar,Primary>& W){
//Update nbLine !!
void UpdateBottomLine(Block<Scalar,Primary>& W1){
int nj = sumNbVectInBlck.back();
//Create temp block
Block<Scalar,Primary> tempRes(nj,W1.getSizeBlock());
call_Tgemm<>(W1.getSizeBlock(),nj,W1.getLeadingDim(),
W1.getPtr(),W1.getLeadingDim(),
getPtrToG(),ldHess,
tempRes.getPtr(),tempRes.getLeadingDim(),
Scalar(1),Scalar(0));
//copy back res inside L (where Wj used to be)
int diff = W1.getLeadingDim()-W1.getSizeBlock();
for(int i=0 ; i<tempRes.getSizeBlock() ; ++i){
memcpy(&getPtrToG()[i*ldHess],&tempRes.getPtr()[i*tempRes.getLeadingDim()],
sizeof(Scalar)*tempRes.getLeadingDim());
//Set to zero the remaing part
memset(&getPtrToG()[i*ldHess + tempRes.getLeadingDim()],0,
sizeof(Scalar)*diff);
}
//Update nbLine
nbLineUsed += W1.getSizeBlock();
}
/**
* @brief Compute Product :
* @brief Compute Least Square and product to get residual :
* |Lj|
* |Hj| x |Y|
*/
template<typename Primary>
void ComputeResidual(Block<Scalar,Primary>& input,
Block<Scalar,Primary>& gamma,
void ComputeResidual(Block<Scalar,Primary>& gamma,
Block<Scalar,Primary>& outputY,
Block<Scalar,Primary>& outputResidual){
//Maybe Gamma need to be augmented
}
......@@ -99,6 +117,7 @@ public:
*
*/
Scalar * getPtrToG(){
return &data[nbLineUsed];
}
/**
......@@ -106,6 +125,8 @@ public:
*
*/
Scalar * getPtrToC(){
//ptrToG + ldHess*n_{j-1}
return getPtrToG() + ldHess*sumNbVectInBlck[sumNbVectInBlck.size()-1];
}
/**
......@@ -113,6 +134,7 @@ public:
*
*/
Scalar * getPtrToD(){
return getPtrToC() + p-nbVectInBlck.back();
}
};
......
......@@ -151,7 +151,7 @@ int main(int ac, char ** av){
Block<Scalar> MatxSolMinusB(SizeBlock,dim);
//X_Diff : X_Exact - Sol_computed
Block<Scalar> X_Diff(SizeBlock,dim);
MatxSolMinusB.CopyBlock(MatxSol,MatxSolMinusB.getSizeBlock());
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];
......
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include <complex>
#include <random>
......@@ -31,5 +32,6 @@ int main(int ac ,char ** av){
}
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