Commit 9c1573cb authored by Cyrille Piacibello's avatar Cyrille Piacibello

Ib should be working soon

parent 53b5d030
......@@ -55,9 +55,8 @@ public:
userEnv);
}
void MatBaseProduct(T * ptrToRead,Block<T,S>& output){
void* toWrite = output.getPtr();
userProduct(userPtrMat,output.getSizeBlock(),ptrToRead,&toWrite,
void MatBaseProduct(T * ptrToRead,int nbRHS,void * ptrToWrite){
userProduct(userPtrMat,nbRHS,ptrToRead,&ptrToWrite,
userEnv);
}
......
......@@ -149,7 +149,7 @@ int main(int ac , char ** av){
double * historic = ibgmresdr_get_logs(handle);
//Write down the convergence
FILE * fd3 = fopen("/home/cyrille/IbGMResDr/matrices_BEM/ConvergenceHisto.txt",
FILE * fd3 = fopen("/home/cyrille/Data/matrices_BEM/ConvergenceHisto.txt",
"w");
fprintf(fd3,"Ite\tMin\tMax\tTime\n");
for(int i=0 ; i<nb_iter_needed ; ++i){
......
......@@ -72,6 +72,15 @@ public:
return currentBlockUsed+1;
}
/**
* @brief return the last block col
*/
Scalar * getLastColPtr(int sizeNewBlock){
Scalar * ptr = getPtr(currentBlockUsed);
addBlock(sizeNewBlock);
return ptr;
}
void addBlockDatas(Block<Scalar,Primary>& inBlock){
//Check
if(inBlock.getLeadingDim() != ldBase){std::cout<<"Problem in addBlockDatas\n";}
......@@ -96,7 +105,10 @@ public:
//Loop over different blocks
for(int i=0 ; i< currentBlockUsed ; ++i){
int sizeCurrBlock = nbVectInBlock[i+1]-nbVectInBlock[i];
std::cout<<"Ortho against block number "<<i<<" of size "
<<sizeCurrBlock<<"\n"
<<"param to gemm : \t"<<sizeCurrBlock<<"\t"<<Wj.getSizeBlock()
<<"\t"<<ldBase<<"\n";
//Compute hessenberg coeff
call_Tgemm<Scalar>(sizeCurrBlock,Wj.getSizeBlock(), ldBase,
getPtr(i),ldBase,
......
......@@ -84,36 +84,52 @@ public:
//QR facto
void ComputeQRFacto(Block& Q, Block& R, //input are already allocated
int vStart=0, int hStart=0){ //if QR of part of block
int hStart=0, int vStart=0){ //if QR of part of block
//Sizes needed
int sizeQrBlock = SizeBlock - hStart;
int sizeQrBlock = SizeBlock - vStart;
//Create Tau
Scalar * tau = new Scalar[sizeQrBlock];
memset(tau,0,sizeof(Scalar)*sizeQrBlock);
//Qr facto can be done in place
std::cout<<"Parameters to call_geqrf : "
<<"\nldb-hStart:\t"<<ldb-hStart
<<"\nsizeQrBlock:\t"<<sizeQrBlock
<<"\nStarting point:\t"<<hStart+vStart*ldb
<<std::endl;
int geqrf = callLapack_geqrf(ldb-hStart,sizeQrBlock,
&datas[vStart+hStart*ldb],ldb,tau);
&datas[hStart+vStart*ldb],ldb,tau);
if(geqrf != 0){
std::cout<<"Problem on QR facto\n";
std::cout<<"Problem on QR facto\n\t output="<<geqrf<<"\n";
}
//Copy triangular R part inside input R
for(int i=0; i<sizeQrBlock ; ++i){
memcpy(R.getPtr(i),&datas[(hStart+i)*ldb+vStart],(i+1)*sizeof(Scalar));
memcpy(R.getPtr(i),&datas[(vStart+i)*ldb+hStart],(i+1)*sizeof(Scalar));
}
//Call to orgqr //orgqr does not exist in Lapack for complex,
//but Lapack has a complex counter part called ungqr, which is
//called if Scalar is a complex type
int orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[vStart+hStart*ldb],ldb,tau);
//Test if caller create a Square space for Q, if so, then we
//use the full QR
int orgqr;
if(Q.getSizeBlock() == Q.getLeadingDim()){
std::cout<<"Generating full Q after QR\n";
orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[hStart+vStart*ldb],ldb,tau,true);
}else{ //Else, only leading p directions are computed
std::cout<<"Generating partial Q after QR\n";
orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[hStart+vStart*ldb],ldb,tau,false);
}
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
}
//Copy output to Q
for(int i=0 ; i<sizeQrBlock ; ++i){
memcpy(Q.getPtr(i),&datas[vStart+(hStart+i)*ldb],
(ldb-vStart)*sizeof(Scalar));
for(int i=0 ; i<Q.getSizeBlock() ; ++i){
memcpy(Q.getPtr(i),&datas[hStart+(vStart+i)*ldb],
(ldb-hStart)*sizeof(Scalar));
}
delete [] tau;
}
......@@ -204,9 +220,14 @@ public:
//Count the number of Singular Value superior to the threshold
//epsilon
int id = 0;
while(sigma[id]>epsilon){
while(id<getSizeBlock() && sigma[id]>epsilon){
id++;
}
//Display all the singular values
for(int i=0 ; i<getSizeBlock() ; ++i){
std::cout<<"sigma["<<i<<"]\t : "<<sigma[i]<<"\n";
}
std::cout<<"Number of sv sup. to "<<epsilon<<"\t: "<<id<<"\n";
output.initData(id,getLeadingDim());
memcpy(output.getPtr(),U.getPtr(),sizeof(Scalar)*id*U.getLeadingDim());
delete [] sigma;
......
#ifndef BLOCKARNOLDIIB_HPP
#define BLOCKARNOLDIIB_HPP
#include "BlockWP.hpp"
#include "HessExtended.hpp"
/**
* @brief : Arnoldi iteration with inexact breakdowns
*
......@@ -12,15 +15,15 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
Primary epsilon = 0.00001){
int SizeBlock = B.getSizeBlock();
int dim = A.size();
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< p ; ++i){
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.getNorm(i));
}
//Compute epsilon_R (step 1)
Primary epsilon_R = normB.front();
for(auto &a : nomrB){
for(auto &a : normB){
if(a<epsilon_R){
epsilon_R = a;
}
......@@ -29,7 +32,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//Here will the restart loop
Base<Scalar> base(p,max_iter+1,dim);
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());
......@@ -63,49 +66,65 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//Main Loop
for(int j=1 ; j<max_iter ; ++j){
//Compute block Wj
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP);
//Compute block Wj inside WP
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
WP.getW());
//Ortho against V and store coeff inside new block column of Hess L
Scalar * toWrite = L.getPtr(j);
base.ComputeMGSOrtho(Wj,toWrite,L.getLeadingDim());
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
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
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);
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
L.ComputeResidual(R1,Y,Res);
std::cout<<"Residual Computed !"<<std::endl;
//Do we really need that ?
std::pair<Primary,Primary> MinMax = getMinMaxNorm();
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm();
//Block to Store U
Block<Scalar,Primary> U;
Res.DecompositionSVD(U,epsilon_R);
std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
<<U.getSizeBlock()<<" !"<<std::endl;
//Create Block to store W
Block<Scalar,Primary> Directions(U.getSizeBlock(),SizeBlock);
Block<Scalar,Primary> tempRpart; //TO BE DISCARDED
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
//TO BE DISCARDED
Block<Scalar,Primary> tempRpart(U.getSizeBlock(),U.getSizeBlock());
U.displayBlock("U: Kept directions");
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
Directions.displayBlock("Q part of U after QR");
//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_{j+1,} = Directions_1^{H} * H_j
//Gj = Directions_2^{H} * H_j
//L_{j+1,} = | Directions_1^{H} | * H_j
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions,U.getSizeBlock());
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;
}
#endif //BLOCKARNOLDIIB_HPP
......@@ -76,7 +76,8 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,Wj);
}else{//Wj = A*V[j] without preCond
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),
Wj.getSizeBlock(),Wj.getPtr());
}
//Where to write in hess
......
......@@ -79,7 +79,11 @@ public:
//3 - Call orgqr on W block
int orgqr = callLapack_orgqr(Parent::getLeadingDim(),getSizeW(),
getW(),Parent::getLeadingDim(),tau);
//check
if(orgqr != 0){
std::cout<<"Problem in BlockWP::QRofW\n";
exit(0);
}
delete [] tau;
}
......@@ -90,8 +94,10 @@ public:
*/
void ComputeProduct(Block<Scalar,Primary>& input,
int nbKeptDir, //Correspond to the number of
//kept directions
Scalar toWrite,int ldToWrite){
//kept directions
Scalar* toWrite,int ldToWrite){
std::cout<<"[P,W] * Dir1 :: nbKeptDir/size : "<<nbKeptDir<<"/"
<<input.getSizeBlock()<<"\n";
call_gemm<>(Parent::getLeadingDim(),nbKeptDir,Parent::getSizeBlock(),
Parent::getPtr(),Parent::getLeadingDim(),
input.getPtr(),input.getLeadingDim(),
......@@ -104,8 +110,8 @@ public:
* write the result in place of P. (obviously used to update P)
*/
void ComputeProduct(Block<Scalar,Primary>& input,
int nbKeptDir){ //Correspond to the start
//of discarded directions
int nbKeptDir){ //Correspond to the start of
//discarded directions
int sizeInputBlock = input.getSizeBlock() - nbKeptDir;
//Create temporary blck to store result,
Block<Scalar,Primary> temp(sizeInputBlock,Parent::getLeadingDim());
......
......@@ -42,12 +42,22 @@ public:
nbLineUsed(0),data(nullptr){
data = new Scalar[ldHess*nbTotalCol];
memset(data,0,sizeof(Scalar)*ldHess*nbTotalCol);
sumNbVectInBlck.push_back(0);
std::cout<<"Init Hess Extended with :\n"
<<"\ttotalSize "<<nbTotalCol*ldHess<<"\n"
<<"\tldHess "<<ldHess<<"\n"
<<"\tSize First block "<<p<<"\n"
<<"\tnbCol "<<nbTotalCol<<"\n";
}
~HessExtended(){
delete [] data;
}
Scalar * getPtr(int idx){
return data;
}
int getLeadingDim() const {
return ldHess;
}
......@@ -56,6 +66,17 @@ public:
return nbLineUsed;
}
void displayHessExtended(std::string name = ""){
std::cout<<name<<"\n";
for(int i=0 ; i<sumNbVectInBlck.back() ; ++i){
for(int j=0 ; j<nbLineUsed+p ; ++j){
std::cout<<data[i*ldHess+j]<<"\t";
}
std::cout<<std::endl;
}
std::cout<<std::endl;
}
/**
* @brief return the ptr to the start of new blockcol, and augment
* the cursors
......@@ -65,36 +86,37 @@ public:
nbVectInBlck.push_back(inSizeBlock);
sumNbVectInBlck.push_back(ref+inSizeBlock);
nbBlckColUsed++;
nbLineUsed += nbVectInBlck.back();
std::cout<<"New Starting col is in data["<<ref*ldHess<<"]\n";
return &data[ref*ldHess];
}
/**
* @brief Remove the last part (containing H), and filling it with
* W^{H}*H
* @brief Update last row of L
*
* @formulae
* * | H_j |
* | L_{j+1,} | | W1^{H} |
* | G_j | = | W2^{H} |
*
* @param W1 Block
* @param W FullBlock
* @param sizeToRead : limite between W1 and W2
*/
template<typename Primary>
void UpdateBottomLine(Block<Scalar,Primary>& W1,int sizeToRead){
void UpdateBottomLine(Block<Scalar,Primary>& W,int sizeToRead){
int nj = sumNbVectInBlck.back();
//Create temp block
Block<Scalar,Primary> tempRes(nj,W1.getSizeBlock());
call_Tgemm<>(sizeToRead,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);
//Save H in order to write directly inside L
Block<Scalar,Primary> Hj(nj,p);
for(int i=0 ; i<nj ; i++){ //Copy H part inside Block Hj
memcpy(Hj.getPtr(i),&getPtrToG()[i*ldHess],
sizeof(Scalar)* p);
}
//Update nbLine
nbLineUsed += sizeToRead;
call_Tgemm(p,nj,p,
W.getPtr(),W.getLeadingDim(),
Hj.getPtr(),Hj.getLeadingDim(),
getPtrToG(),ldHess,
Scalar(1),Scalar(0));
//nbLineUsed += sizeToRead;
}
/**
......@@ -107,13 +129,20 @@ public:
Block<Scalar,Primary>& outputY,
Block<Scalar,Primary>& outputResidual){
//Complete gamma
std::cout<<"About to compute Residual...\n";
std::cout<<"Gamma : is "<<gamma.getSizeBlock()<<" x "
<<gamma.getLeadingDim()<<"\n";
std::cout<<"Nb Line used : "<<nbLineUsed<<"\n";
displayHessExtended("Hess before gels");
int ldR = nbLineUsed + p;
Scalar * r = new Scalar[ldR*p];
memset(r,0,((nbLineUsed + p)*p)*sizeof(Scalar));
std::cout<<"R : is "<<p<<" x "<<ldR<<"...\n";
for(int i=0 ; i<p ; ++i){
memcpy(&r[i*ldR],&gamma.getPtr(i)[i*gamma.getLeadingDim()],
memcpy(&r[i*ldR],gamma.getPtr(i),
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));
......@@ -125,6 +154,9 @@ public:
int gels = callLapack_gels<>(ldR,sumNbVectInBlck.back(),p,
workingCopy,ldR,
r,ldR);
if(gels != 0){
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
......@@ -132,13 +164,16 @@ public:
data,ldHess,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(-1),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.displayBlock("Residual obtained from LeastSquare");
delete [] workingCopy;
delete [] r;
}
......@@ -153,6 +188,7 @@ public:
*
*/
Scalar * getPtrToG(){
// std::cout<<"\tCall to getPtrToG():\n\t return data["<<nbLineUsed<<"]\n";
return &data[nbLineUsed];
}
......@@ -162,7 +198,9 @@ public:
*/
Scalar * getPtrToC(){
//ptrToG + ldHess*n_{j-1}
return getPtrToG() + ldHess*sumNbVectInBlck[sumNbVectInBlck.size()-1];
// std::cout<<"\tCall to getPtrToC():\n\t\t return getPtrToG()+"
// <<ldHess<<" * "<<sumNbVectInBlck[sumNbVectInBlck.size()-2]<<"\n";
return getPtrToG() + ldHess*sumNbVectInBlck[sumNbVectInBlck.size()-2];
}
/**
......@@ -170,6 +208,8 @@ public:
*
*/
Scalar * getPtrToD(){
// std::cout<<"\tCall to getPtrToD():\n\t\t\t return getPtrToC()+"
// <<p<<" - "<<nbVectInBlck.back()<<"\n";
return getPtrToC() + p-nbVectInBlck.back();
}
......
......@@ -238,6 +238,11 @@ int callLapack_gels<std::complex<double> >(int nbMaxEl, int nbCol, int nbRHS,
/**
* QR-Factorisation (full)
* @param m number of lines
* @param n number of col
* @param a : data to factorize
* @param lda : leading dim of a
* @param tau : array of size at least min(m,n) (output)
*/
template<class Scalar>
int callLapack_geqrf(int m,int n,Scalar* a,int lda,Scalar * tau){
......@@ -276,33 +281,57 @@ int callLapack_geqrf<std::complex<double> >(int m,int n,
* The only Q part generated is the orthogonal part (spanning of
* Im(A))
*
* @param m Order of matrix Q ?
* @param p Nb vector needed
*
* doc : https://software.intel.com/en-us/node/521010
*/
template<typename Scalar>
int callLapack_orgqr(int , int , Scalar * ,int , Scalar *){
int callLapack_orgqr(int , int , Scalar * ,int , Scalar *,bool=false){
std::cout<<"Only Float or Double supported (orgqr)"<< std::endl;
exit(0);
return 1;
}
template<>
int callLapack_orgqr(int m, int p, float * mat, int ldMat, float *tau){
return LAPACKE_sorgqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
int callLapack_orgqr(int m, int p,float * mat, int ldMat, float *tau,
bool full){
if(full){
return LAPACKE_sorgqr(LAPACK_COL_MAJOR, m, m, p, mat, ldMat,tau);
}else{
return LAPACKE_sorgqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
}
}
template<>
int callLapack_orgqr(int m, int p, double * mat, int ldMat, double *tau){
return LAPACKE_dorgqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
int callLapack_orgqr(int m, int p, double * mat, int ldMat, double *tau,
bool full){
if(full){
return LAPACKE_dorgqr(LAPACK_COL_MAJOR, m, m, p, mat, ldMat,tau);
}else{
return LAPACKE_dorgqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
}
}
//https://software.intel.com/en-us/node/521012#FCFCCE90-1C8F-4879-824A-7EBA7C119AE8
template<>
int callLapack_orgqr(int m, int p,
std::complex<float> * mat, int ldMat,
std::complex<float> * tau){
return LAPACKE_cungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
std::complex<float> * tau, bool full){
if(full){
return LAPACKE_cungqr(LAPACK_COL_MAJOR, m, m, p, mat, ldMat,tau);
}else{
return LAPACKE_cungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
}
}
template<>
int callLapack_orgqr(int m, int p,
std::complex<double> * mat, int ldMat,
std::complex<double> * tau){
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
std::complex<double> * tau, bool full){
if(full){
std::cout<<"There !\n"<<
"Dans l'ordre : m "<<m<<" p "<<p<<"\n";;
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, m, p, mat, ldMat,tau);
}else{
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
}
}
/**
......
......@@ -6,6 +6,7 @@ set(GMRES_TEST_SRC
testBlockGMRes.cpp
testBlockSvd.cpp
testBlockWP.cpp
testBlockArnoldiIb.cpp
)
#set(LIBS_FOR_TESTS ibgmresdr)
......
......@@ -56,11 +56,11 @@ public:
Scalar(1),Scalar(0));
}
void MatBaseProduct(Scalar * toRead,Block<Scalar>& out){
call_gemm<Scalar>(dim,out.getSizeBlock(), dim,
void MatBaseProduct(Scalar * toRead,int nbRHS, Scalar * ptrToWrite){
call_gemm<Scalar>(dim,nbRHS, dim,
data,dim,
toRead,dim,
out.getPtr(),out.getLeadingDim(),
ptrToWrite,dim,
Scalar(1),Scalar(0));
}
......@@ -130,7 +130,7 @@ int main(int ac, char ** av){
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
X_Exact.getPtr(i)[j] = dis(gen);
X0.getPtr(i)[j] = dis(gen);
}
}
......
......@@ -22,7 +22,7 @@ int main(int ac, char ** av){
//Fill the block
for(int i=0; i<SizeBlock ; ++i){
for(int j=0 ; j<sizeVect ; ++j){
B.getPtr()[ i* B.getLeadingDim() + j] = dis(gen);
B.getPtr()[ i* B.getLeadingDim() + j] = {dis(gen),dis(gen)};
}
}
......@@ -34,6 +34,8 @@ int main(int ac, char ** av){
B.DecompositionSVD(U,epsilon);
U.displayBlock("Left singular vectors : U");
std::cout<<"Number of singular values kept over max number : "
<<U.getSizeBlock()<<"/"<<SizeBlock<<std::endl;
......
......@@ -44,7 +44,7 @@ int main(int ac ,char ** av){
for(int j=0 ; j<Dim ; ++j){ //Fill it
v.getPtr()[j] = {dis(gen),dis(gen)};
}
WP.ComputeProduct(v);
WP.ComputeProduct(v,0);
WP.displayBlock();
Block<Scalar,Primary> C(WP.getSizeW(),WP.getSizeP());
......
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