Commit 96b8debc authored by Cyrille Piacibello's avatar Cyrille Piacibello

testBlock updated and working testHessExtended to be added

parent b3976813
......@@ -28,10 +28,10 @@ public:
return Parent::getPtr();
}
int getSizeP() const {
int getSizeP(){
return cursor;
}
int getSizeW() const {
int getSizeW(){
return Parent::getSizeBlock()-cursor;
}
......@@ -41,7 +41,7 @@ public:
*/
void OrthoWP(Scalar * toWrite, int ldToWrite){
if(cursor>0){
//Generate coeef inside C
//Generate coef inside C
call_Tgemm<>(getSizeP(),getSizeW(),Parent::getLeadingDim(),
getP(),Parent::getLeadingDim(),
getW(),Parent::getLeadingDim(),
......@@ -114,8 +114,10 @@ public:
//updated
memcpy(getP(),temp.getPtr(),
sizeof(Scalar)*temp.getSizeBlock()*temp.getLeadingDim());
cursor = temp.getSizeBlock();
//Not necessary, should be erased but still...
//Not necessary, should be written over at the next iteration
//but still...
memset(getW(),0,getSizeW()*Parent::getLeadingDim()*sizeof(Scalar));
}
};
......
......@@ -102,7 +102,39 @@ public:
void ComputeResidual(Block<Scalar,Primary>& gamma,
Block<Scalar,Primary>& outputY,
Block<Scalar,Primary>& outputResidual){
//Complete gamma
int ldR = nbLineUsed + p;
Scalar * r = new Scalar[ldR*p];
memset(r,0,((nbLineUsed + p)*p)*sizeof(Scalar));
for(int i=0 ; i<p ; ++i){
memcpy(&r[i*ldR],&gamma.getPtr(i)[i*gamma.getLeadingDim()],
sizeof(Scalar)*(i+1));
}
//get a working copy of F
Scalar*workingCopy = new Scalar[ldR*sumNbVectInBlck.back()];
memset(workingCopy,0,ldR*sumNbVectInBlck.back()*sizeof(Scalar));
for(int i=0 ; i<sumNbVectInBlck.back() ; ++i){
memcpy(&workingCopy[i*ldR],&data[i*ldHess],
sizeof(Scalar)*ldR);
}
//Call least square
int gels = callLapack_gels<>(ldR,sumNbVectInBlck.back(),p,
workingCopy,ldR,
r,ldR);
//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(),
data,ldHess,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(-1),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];
}
}
}
......
......@@ -186,7 +186,7 @@ int callLapack_gels<double>(int nbMaxEl, int nbCol, int nbRHS,
return LAPACKE_dgels(LAPACK_COL_MAJOR, //column major
'N', //Something {Nope,Trans,Hermittienne}
nbMaxEl, //number of elements inside column
nbCol, //number of elements in each columns
nbCol, //number of columns
nbRHS, //number of rhs
matrice, //Matrice
ldMatrice, //leading dimension
......
......@@ -12,8 +12,8 @@ int main(int ac ,char ** av){
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock = 5;
int Dim = 25;
int SizeBlock = 3;
int Dim = 10;
using Primary=double;
using Scalar=std::complex<Primary>;
......@@ -24,14 +24,40 @@ int main(int ac ,char ** av){
BlockWP<Scalar,Primary> WP(SizeBlock,Dim);
//Fill it
// Fill it
for(int i=0; i<SizeBlock ; ++i){
for(int j=0 ; j<Dim ; ++j){
WP.getPtr()[ i* WP.getLeadingDim() + j] = dis(gen);
WP.getPtr(i)[j] = {dis(gen),dis(gen)};
}
}
WP.displayBlock("Initial WP");
//Create block
Block<Scalar,Primary> R(3,3);
WP.QRofW(R.getPtr(),R.getLeadingDim());
WP.displayBlock("WP after QR");
R.displayBlock("R part");
//Create vector
Block<Scalar,Primary> v(1,Dim);
for(int j=0 ; j<Dim ; ++j){ //Fill it
v.getPtr()[j] = {dis(gen),dis(gen)};
}
WP.ComputeProduct(v);
WP.displayBlock();
Block<Scalar,Primary> C(WP.getSizeW(),WP.getSizeP());
WP.OrthoWP(C.getPtr(),C.getLeadingDim());
//Check ortho
for(int i = 0 ; i<WP.getSizeW() ; ++i){
Scalar acc = 0;
for(int j=0 ; j<WP.getLeadingDim() ; ++j){
acc += WP.getP()[j] * WP.getW()[i*WP.getLeadingDim() + j];
}
std::cout<<"Ortho between P and W["<<i<<"] : "<<acc<<"\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