Commit 9bf2d49b authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

remove MORSE_enable(MORSE_ERROR) (not supported anymore in recent version of chameleon)

parent c48f0dfa
......@@ -94,7 +94,6 @@ static double* create_random_matrix(int M, int N)
int main(int argc, char *argv[])
{
MORSE_Init(2, 0);
MORSE_Enable(MORSE_ERRORS);
int dim = 25;
if (argc > 1)
......
......@@ -166,7 +166,7 @@ struct Arnoldi_IB
Block LS; // Block to hold LS
// STEP1: Compute first block residual
Block R0{SizeBlock,dim};
Block R0{SizeBlock, dim};
//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
......@@ -220,6 +220,7 @@ struct Arnoldi_IB
// Ortho and filling of L (the (IB) Hessenberg)
Arnoldi_Ortho(L, base, WP, ArnChoice, ortho, A);
//L.displayHessExtendedBitMap("AFTER ORTHO");
// Create Block and Residual
LS.initData(SizeBlock, L.getNbLineUsed() + SizeBlock);
......@@ -292,6 +293,7 @@ struct Arnoldi_IB
std::cout<<"################## End of Iteration ##################\n\n";
log.add_ite(ite_meter,WP.getSizeW(),evalConv.Min,evalConv.Max,timestamp);
j += WP.getSizeW();
//L.displayHessExtendedBitMap("END L HESS");
}
if (!convergence) { //No convergence, we need to write current sol into X0
Block Sol{SizeBlock,dim};
......
......@@ -16,13 +16,13 @@ void ArnoldiBlock(Hess &L, Base &base, BlockWP &WP, OrthoChoice ortho)
// We already test if there was room to augment Hessenberg
auto *ptrToHess = L.getNewStartingCol(WP.getSizeW());
//Ortho against vectors inside base
// Ortho against vectors inside base
base.ComputeOrtho(ortho, WP, ptrToHess, L.getLeadingDim(), WP.getSizeP());
std::cout << "Ortho done" << std::endl;
// Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(), L.getLeadingDim());
WP.QRofW(L.getPtrToD(),L.getLeadingDim()); // QR of w part
WP.QRofW(L.getPtrToD(), L.getLeadingDim()); // QR of w part
/* Warning : No base augmentation because we need to know if IB
* already happened at this stage */
......
......@@ -23,22 +23,20 @@ void ArnoldiRuhe_CGS(Hess& H,
BlockWP<Scalar,Primary,KI>& WP,
Matrix &A )
{
//Cste
int dim = base.getLeadingDim();
//Where to start writing into Hessenberg
const int dim = base.getLeadingDim();
// Where to start writing into Hessenberg
Scalar *ptrToHess = H.getNewStartingCol(WP.getSizeW());
//Loop over vector in WP block
for (int k=0; k< WP.getSizeW(); ++k) {
int nbVectInBase = base.getNbVectUsed();
//Starting ptr of W_k
// Starting ptr of W_k
Scalar *ptrToWK = WP.getW() + k * WP.getLeadingDim();
//Starting point of Hessenberg
// Starting point of Hessenberg
Scalar *ptrToHessK = ptrToHess + k * H.getLeadingDim();
//Starting point of C inside Hessenberg
// Starting point of C inside Hessenberg
Scalar *ptrToCK = H.getPtrToC() + k * H.getLeadingDim();
//Starting point of D inside Hessenberg
// Starting point of D inside Hessenberg
Scalar *ptrToDK = H.getPtrToD() + k * H.getLeadingDim();
{//Compute Ortho against base
......
......@@ -57,12 +57,12 @@ int BGMRes(Matrix &A,
print_precond(A);
std::vector<Primary> normB, invNormB;
B.Normalize(normB,A); // Norm(s) of B
B.Normalize(normB, A); // Norm(s) of B
std::transform(normB.begin(),normB.end(),std::back_inserter(invNormB),
[](const Primary &a)-> Primary{ return 1/a; });
X0.ScaleBlock(invNormB);
// Global step is total number of iteration (thus is <max_iter)
// Global step is total number of iteration (thus is < max_iter)
int nb_restart = 0, nbIteDone = 0, KrylovSpaceSpanned = 0;
ArnReturn ret{0, 0, false};
......@@ -79,7 +79,6 @@ int BGMRes(Matrix &A,
Base base{B.getSizeBlock(),SizeToSpan,dim};
Block LS;
DRParam<Scalar,Primary> DR_struct{false,hessSave,base,0,Scalar(0),LS};
ret = ARNOLDI::Exec_Procedure( // Call Arnoldi's procedure specified in template
A, X0, B, SizeToSpan, epsilon,log,DR_struct, ortho, arnChoice);
......
......@@ -82,10 +82,8 @@ public:
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*ld],
(i+1)*sizeof(Scalar));
}
for (int i=0; i<getSizeW(); ++i)
memcpy(&toWriteR[i*ldR], &getW()[i*ld], (i+1)*sizeof(Scalar));
//3 - Call orgqr on W block
Block<Scalar,Primary,KI> tmp{getSizeW(), ld};
......
......@@ -193,10 +193,10 @@ public:
sizeof(Scalar)* _SizeBlock);
}
KI::Tgemm(_SizeBlock,nj,_SizeBlock,
W.getPtr(),W.getLeadingDim(),
Hj.getPtr(),Hj.getLeadingDim(),
getPtrToG(),_ld,
Scalar(1),Scalar(0));
W.getPtr(),W.getLeadingDim(),
Hj.getPtr(),Hj.getLeadingDim(),
getPtrToG(),_ld,
Scalar(1),Scalar(0));
}
/**
......@@ -216,7 +216,7 @@ public:
// Create New Phi
Block<Scalar,Primary,KI> * newPhi =
new Block<Scalar,Primary,KI>{_SizeBlock,
Phi->getLeadingDim() + p_jplus1};
Phi->getLeadingDim() + p_jplus1};
// Init first n_j rows
for (int i=0; i<_SizeBlock; ++i)
......@@ -224,11 +224,11 @@ public:
// Compute W1W2^{H} * Phi(nj+1:nj+p,:)
KI::Tgemm(_SizeBlock,_SizeBlock,
_SizeBlock,
Directions.getPtr(),Directions.getLeadingDim(),
&Phi->getPtr()[nj],Phi->getLeadingDim(),
&newPhi->getPtr()[nj],newPhi->getLeadingDim(),
Scalar(1),Scalar(0));
_SizeBlock,
Directions.getPtr(),Directions.getLeadingDim(),
&Phi->getPtr()[nj],Phi->getLeadingDim(),
&newPhi->getPtr()[nj],newPhi->getLeadingDim(),
Scalar(1),Scalar(0));
delete Phi;
Phi = newPhi;
}
......@@ -241,10 +241,10 @@ public:
void ComputeRHS(Scalar * RHS, int ldRHS)
{
KI::gemm(Phi->getLeadingDim(),_SizeBlock,_SizeBlock,
Phi->getPtr(),Phi->getLeadingDim(),
Lambda1->getPtr(),Lambda1->getLeadingDim(),
RHS,ldRHS,
Scalar(1),Scalar(0));
Phi->getPtr(),Phi->getLeadingDim(),
Lambda1->getPtr(),Lambda1->getLeadingDim(),
RHS,ldRHS,
Scalar(1),Scalar(0));
std::cout<<"Computing RHS done at ite "<< sumNbVectInBlock.size() <<"\n";
}
......@@ -258,55 +258,50 @@ public:
{
// Complete gamma
int ldR = nbLineUsed + _SizeBlock;
Scalar * r = new Scalar[ldR*_SizeBlock];
memset(r,0,ldR*_SizeBlock*sizeof(Scalar));
if (!Phi) {
for (int i=0; i<_SizeBlock; ++i) {
memcpy(&r[i*ldR],Lambda1->getPtr(i),
sizeof(Scalar)*(i+1));
}
} else {
int nj = sumNbVectInBlock.back();
Scalar *r = new Scalar[ldR*_SizeBlock];
Scalar *L_copy = new Scalar[ldR*nj];
memset(r, 0, ldR*_SizeBlock*sizeof(Scalar));
// r = Λ
if (!Phi)
for (int i=0; i<_SizeBlock; ++i)
memcpy(&r[i*ldR],Lambda1->getPtr(i), sizeof(Scalar)*(i+1));
else
ComputeRHS(r, ldR);
}
Block<Scalar,Primary,KI> RHS{
_SizeBlock,outputResidual.getLeadingDim()};
Block<Scalar,Primary,KI> RHS{_SizeBlock, outputResidual.getLeadingDim()};
// RHS = r ( = Λ )
for (int i=0; i<outputResidual.getSizeBlock();++i)
for (int j=0; j<outputResidual.getLeadingDim(); ++j)
RHS.getPtr(i)[j] = r[i*ldR+j];
// get a working copy of F
// displayHessExtendedBitMap("Before calling gels");
// displayHessExtendedBitMap("Before calling gels");
// displayHessExtended("Same");
Scalar *workingCopy = new Scalar[ldR*sumNbVectInBlock.back()];
memset(workingCopy, 0, ldR*sumNbVectInBlock.back()*sizeof(Scalar));
for (int i=0; i<sumNbVectInBlock.back(); ++i) {
memcpy(&workingCopy[i*ldR],
&_data[i*_ld],
sizeof(Scalar)*ldR);
}
memset(L_copy, 0, ldR*nj*sizeof(Scalar));
for (int i=0; i<nj; ++i) // L_copy = L
memcpy(&L_copy[i*ldR], &_data[i*_ld], sizeof(Scalar)*ldR);
// displayHessExtendedBitMap("Residual");
// solve least square
int gels = KI::gels(ldR, sumNbVectInBlock.back(), _SizeBlock,
workingCopy, ldR,
r, ldR);
int gels = KI::gels(ldR, nj, _SizeBlock, L_copy, ldR, r, ldR);
if (gels != 0)
std::cout<<"Problem on Least Square\n\t gels="<<gels<<"\n";
//Copy back res inside Y
for (int i=0; i<_SizeBlock; ++i)
memcpy(outputY.getPtr(i),&r[i*ldR],sizeof(Scalar)*sumNbVectInBlock.back());
memcpy(outputY.getPtr(i),&r[i*ldR], sizeof(Scalar)*nj);
//Call gemm : outputRes = F*Y-gama
KI::gemm(nbLineUsed + _SizeBlock,
_SizeBlock,sumNbVectInBlock.back(),
_data,_ld,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(1),Scalar(0));
KI::gemm(nbLineUsed + _SizeBlock, _SizeBlock, nj,
_data, _ld,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(1),Scalar(0));
std::cout << "Compute F*Y done\n";
......@@ -315,7 +310,7 @@ public:
for (int j=0; j<outputResidual.getLeadingDim(); ++j)
outputResidual.getPtr(i)[j] -= RHS.getPtr(i)[j];
delete [] workingCopy;
delete [] L_copy;
delete [] r;
}
......
......@@ -10,7 +10,6 @@ int main(int argc, char *argv[])
std::cout << arg << "\n";
MORSE_Init(2, 0);
MORSE_Enable(MORSE_ERRORS);
using Primary = double;
using Scalar = std::complex<Primary>;
......@@ -25,7 +24,8 @@ int main(int argc, char *argv[])
mml.LoadMatrix(mat, filename);
int dim = mat.size();
int SizeBlock = 16, restart = 0, nbEigenVector = 6, maxMVP = 2000;;
int SizeBlock = 16, restart = 0, nbEigenVector = 6, maxMVP = 2000;
if (args.size() > 1)
SizeBlock = std::stoi(args[1]);
if (args.size() > 2)
......@@ -38,7 +38,7 @@ int main(int argc, char *argv[])
Block XExact{SizeBlock, dim};
Block RHS{SizeBlock, dim};
std::vector<Primary> epsilon{1e-5};
std::vector<Primary> epsilon{1e-4};
for (int i=0; i<SizeBlock; ++i)
XExact.getPtr(i)[i] = Scalar{1.0/*,1.0*/};
......@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
"young1c",
mat, RHS, XExact,
restart, maxMVP, epsilon,
OrthoChoice::CGS, ArnOrtho::BLOCK
OrthoChoice::IMGS, ArnOrtho::BLOCK
);
// auto r2 = runTest_BGMRES_filelog<Arnoldi_IB>(
// "young1c",
......
......@@ -10,7 +10,6 @@ int main(int argc, char *argv[])
std::cout << arg << "\n";
MORSE_Init(2, 0);
MORSE_Enable(MORSE_ERRORS);
MORSE_user_tag_size(31, 16);
int SizeBlock = 16;
......
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