CDMATH/linearsolver/src/LinearSolver.cxx

   1 /**
   2  * LinearSolver.cxx
   3  *
   4  *  Created on: 15 avr. 2013
   5  *      Author: mekkas
   6  */
   7
   8 #include <cmath>
   9 #include <string>
  10
  11 #include <petscksp.h>
  12 #include <petscmat.h>
  13 #include <petscvec.h>
  14
  15 #include "CdmathException.hxx"
  16 #include "LinearSolver.hxx"
  17 #include "SparseMatrixPetsc.hxx"
  18
  19 using namespace std;
  20
  21
  22 LinearSolver::LinearSolver ( void )
  23 {
  24         _tol=1.E-15;
  25         _numberMaxOfIter=0;
  26         _residu=1.E30;
  27         _convergence=false;
  28         _numberOfIter=0;
  29         _isSingular=false;
  30         _nameOfPc="";
  31         _nameOfMethod="";
  32         _mat=NULL;
  33         _smb=NULL;
  34         _ksp=NULL;
  35         _prec=NULL;
  36         _isSparseMatrix=false;
  37         _computeConditionNumber=false;
  38     _conditionNumber=1e100;
  39 }
  40
  41 void
  42 LinearSolver::kspDuplicate(const KSP source, const Mat mat, KSP* destination) const
  43 {
  44         KSPCreate(PETSC_COMM_WORLD,&(*destination));
  45 #ifdef PETSC_VERSION_GREATER_3_5
  46         KSPSetOperators(*destination,mat,mat);
  47 #else
  48         KSPSetOperators(*destination,mat,mat,SAME_NONZERO_PATTERN);
  49 #endif
  50         KSPType type;
  51         KSPGetType(source,&type);
  52         KSPSetType(*destination,type);
  53         /*
  54     PetscReal tol1,tol2,tol3;
  55     PetscInt maxIter;
  56     KSPGetTolerances(source,&tol1,&tol2,&tol3,&maxIter);
  57     KSPSetTolerances(*destination,tol1,tol2,tol3,maxIter);
  58         // */
  59 }
  60
  61 void
  62 LinearSolver::precDuplicate(const PC source, const KSP ksp, PC* destination)  const
  63 {
  64         KSPGetPC(ksp,destination);
  65         PCType type;
  66         PCGetType(source,&type);
  67         PCSetType(*destination,type);
  68 }
  69
  70 LinearSolver::~LinearSolver ( void )
  71 {
  72         //if(&_mat != NULL)
  73         //      MatDestroy(&_mat);
  74         if(&_smb != NULL)
  75                 VecDestroy(&_smb);
  76         KSPDestroy(&_ksp);
  77         //PetscFinalize();
  78 }
  79
  80 void
  81 LinearSolver::setTolerance(double tol)
  82 {
  83         _tol=tol;
  84         KSPSetTolerances(_ksp,tol,PETSC_DEFAULT,PETSC_DEFAULT,getNumberMaxOfIter());
  85 }
  86
  87 void
  88 LinearSolver::setNumberMaxOfIter(int numberMaxOfIter)
  89 {
  90         _numberMaxOfIter=numberMaxOfIter;
  91         KSPSetTolerances(_ksp,getTolerance(),PETSC_DEFAULT,PETSC_DEFAULT,numberMaxOfIter);
  92 }
  93
  94 void
  95 LinearSolver::setComputeConditionNumber(bool display)
  96 {
  97         _computeConditionNumber=display;
  98 }
  99
 100 double
 101 LinearSolver::getConditionNumber() const
 102 {
 103     if(_computeConditionNumber and _convergence and !_isSingular)
 104         return _conditionNumber;
 105     else if(!_computeConditionNumber)
 106         throw CdmathException("LinearSolver::getConditionNumber(): Condition number can not be evaluated without prior call to setComputeConditionNumber()");
 107     else if(_isSingular)//matrix is singular
 108         throw CdmathException("LinearSolver::getConditionNumber(): Condition number can not be evaluated by PETSC for singular matrices (division by zero)");
 109     else //_convergence = false
 110         throw CdmathException("LinearSolver::getConditionNumber(): Condition number can not be evaluated without prior successful resolution of a linear system");
 111 }
 112
 113 LinearSolver::LinearSolver( const GenericMatrix& matrix,
 114                 const Vector& secondMember,
 115                 int numberMaxOfIter,
 116                 double tol,
 117                 string nameOfMethod,
 118                 string nameOfPc )
 119 {
 120         _tol = tol;
 121         _numberMaxOfIter = numberMaxOfIter;
 122         _residu = 1.E30;
 123         _convergence = false;
 124         _numberOfIter = 0;
 125         _isSingular = false;
 126         _isSparseMatrix = matrix.isSparseMatrix();
 127         _computeConditionNumber=false;
 128         _nameOfPc = nameOfPc;
 129         _nameOfMethod = nameOfMethod;
 130         _secondMember = secondMember;
 131         _mat = NULL;
 132         _smb = NULL;
 133         _prec = NULL;
 134         _ksp = NULL;
 135
 136         //setTolerance(tol);
 137         //setNumberMaxOfIter(numberMaxOfIter);
 138         setPreconditioner(nameOfPc);
 139         setMethod(nameOfMethod);
 140         setLinearSolver(matrix, secondMember);
 141 }
 142
 143 LinearSolver::LinearSolver( const GenericMatrix& matrix,
 144                 const std::vector<double>& secondMember,
 145                 int numberMaxOfIter,
 146                 double tol,
 147                 string nameOfMethod,
 148                 string nameOfPc )
 149 {
 150         _tol = tol;
 151         _numberMaxOfIter = numberMaxOfIter;
 152         _residu = 1.E30;
 153         _convergence = false;
 154         _numberOfIter = 0;
 155         _isSingular = false;
 156         _isSparseMatrix = matrix.isSparseMatrix();
 157         _computeConditionNumber=false;
 158         _nameOfPc = nameOfPc;
 159         _nameOfMethod = nameOfMethod;
 160         _mat = NULL;
 161         _smb = NULL;
 162         _prec = NULL;
 163         _ksp = NULL;
 164
 165         //setTolerance(tol);
 166         //setNumberMaxOfIter(numberMaxOfIter);
 167         setPreconditioner(nameOfPc);
 168         setMethod(nameOfMethod);
 169     //converting vector<double> to Vector
 170     int size=secondMember.size();
 171     Vector secondMemberVector(size);
 172         for ( int i=0; i<size; i++)
 173                 secondMemberVector(i)=secondMember[i];
 174         _secondMember = secondMemberVector;
 175
 176         setLinearSolver(matrix, secondMemberVector);
 177 }
 178
 179 void
 180 LinearSolver::setPreconditioner(string pc)
 181 {
 182         if((pc.compare("ICC") != 0) && (pc.compare("ILU") != 0) && (pc.compare("LU") != 0) && (pc.compare("CHOLESKY") != 0) && (pc.compare("") != 0))
 183         {
 184                 string msg="LinearSolver::LinearSolver: preconditioner "+pc+" does not exist.\n";
 185                 throw CdmathException(msg);
 186         }
 187         if ((pc.compare("ICC")==0 || pc.compare("ILU")==0) && _isSparseMatrix==false)
 188         {
 189                 string msg="LinearSolver::LinearSolver: preconditioner "+pc+" is not compatible with dense matrix.\n";
 190                 throw CdmathException(msg);
 191         }
 192
 193         if ((pc.compare("ICC")==0 || pc.compare("ILU")==0) && (_nameOfMethod.compare("LU")==0 || _nameOfMethod.compare("CHOLESKY")==0 ))
 194         {
 195                 string msg="LinearSolver::LinearSolver: preconditioner "+pc+" is not compatible with "+_nameOfMethod+".\n";
 196                 throw CdmathException(msg);
 197         }
 198         _nameOfPc=pc;
 199 }
 200
 201
 202 void
 203 LinearSolver::setMethod(string nameOfMethod)
 204 {
 205         _nameOfMethod = nameOfMethod;
 206
 207         if ((_nameOfPc.compare("ICC")==0 || _nameOfPc.compare("ILU")==0) && (_nameOfMethod.compare("LU")==0 || _nameOfMethod.compare("CHOLESKY")==0) )
 208         {
 209                 string msg="LinearSolver::LinearSolver: preconditioner "+_nameOfPc+" is not compatible with "+_nameOfMethod+".\n";
 210                 throw CdmathException(msg);
 211         }
 212 }
 213
 214
 215 void
 216 LinearSolver::setLinearSolver(const GenericMatrix& matrix, const Vector& secondMember)
 217 {
 218         if ((_nameOfPc.compare("ICC")==0 || _nameOfPc.compare("ILU")==0) && _isSparseMatrix==false)
 219         {
 220                 string msg="LinearSolver::LinearSolver: preconditioner "+_nameOfPc+" is not compatible with dense matrix.\n";
 221                 throw CdmathException(msg);
 222         }
 223
 224         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);//All constructors lead here so we initialize petsc here
 225         setMatrix(matrix);
 226         setSndMember(secondMember);
 227 }
 228
 229
 230 bool
 231 LinearSolver::isSparseMatrix( void ) const
 232 {
 233         return (_isSparseMatrix);
 234 }
 235
 236 bool
 237 LinearSolver::isMatrixSingular( void ) const
 238 {
 239         return _isSingular;
 240 }
 241
 242 int
 243 LinearSolver::getNumberOfIter( void ) const
 244 {
 245         return (_numberOfIter);
 246 }
 247
 248 bool
 249 LinearSolver::getStatus( void ) const
 250 {
 251         return (_convergence);
 252 }
 253
 254 double
 255 LinearSolver::getResidu( void ) const
 256 {
 257         return (_residu);
 258 }
 259
 260 double
 261 LinearSolver::getTolerance(void) const
 262 {
 263         return (_tol);
 264 }
 265
 266 int
 267 LinearSolver::getNumberMaxOfIter(void) const
 268 {
 269         return (_numberMaxOfIter);
 270 }
 271
 272 void
 273 LinearSolver::setMatrix(const GenericMatrix& matrix)
 274 {
 275         /* matrix to mat */
 276         int numberOfRows=matrix.getNumberOfRows();
 277         int numberOfColumns=matrix.getNumberOfColumns();
 278
 279         if (matrix.isSparseMatrix())//Matrix is already a petsc structure
 280         {
 281                         const SparseMatrixPetsc& Smat = dynamic_cast<const SparseMatrixPetsc&>(matrix);
 282                         _mat=Smat.getPetscMatrix();
 283         }
 284         else//Matrix is in A. Mekkas dense structure
 285         {
 286                 MatCreate(PETSC_COMM_WORLD, &_mat);
 287                 MatSetSizes(_mat, PETSC_DECIDE, PETSC_DECIDE, numberOfRows, numberOfColumns);
 288                 MatSetType(_mat,MATSEQDENSE);
 289
 290                 PetscScalar *a;
 291                 PetscMalloc(numberOfRows*numberOfColumns*sizeof(PetscScalar),&a);
 292                 MatSeqDenseSetPreallocation(_mat,a);
 293
 294                 for (int i=0;i<numberOfRows;i++)
 295                         for (int j=0;j<numberOfColumns;j++)
 296                                 a[i+j*numberOfRows]=matrix(i,j);
 297         }
 298         //Assemblage final
 299         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 300         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 301
 302         KSPCreate(PETSC_COMM_WORLD, &_ksp);
 303 #ifdef PETSC_VERSION_GREATER_3_5
 304 KSPSetOperators(_ksp,_mat,_mat);
 305 #else
 306         KSPSetOperators(_ksp,_mat,_mat,SAME_NONZERO_PATTERN);
 307 #endif
 308
 309         KSPGetPC(_ksp,&_prec);
 310 }
 311
 312 void
 313 LinearSolver::setSndMember(const Vector& secondMember)
 314 {
 315         _secondMember=secondMember;
 316         if(&_smb!=NULL)
 317                 VecDestroy(&_smb);
 318         _smb=vectorToVec(secondMember);
 319
 320 }
 321
 322 void
 323 LinearSolver::setMatrixIsSingular(bool sing)
 324 {
 325         _isSingular=sing;
 326 }
 327
 328 Vector
 329 LinearSolver::getSndMember(void) const
 330 {
 331         return (_secondMember);
 332 }
 333
 334 string
 335 LinearSolver::getNameOfMethod(void) const
 336 {
 337         return (_nameOfMethod);
 338 }
 339
 340 string
 341 LinearSolver::getNameOfPc(void) const
 342 {
 343         return (_nameOfPc);
 344 }
 345
 346 LinearSolver::LinearSolver ( const LinearSolver& LS )
 347 {
 348         _tol=LS.getTolerance();
 349         _nameOfMethod=LS.getNameOfMethod();
 350         _numberMaxOfIter=LS.getNumberMaxOfIter();
 351         _secondMember=LS.getSndMember();
 352         _residu=LS.getResidu();
 353         _convergence=LS.getStatus();
 354         _numberOfIter=LS.getNumberOfIter();
 355         _isSingular=LS.isMatrixSingular();
 356         _nameOfPc=LS.getNameOfPc();
 357         _mat=NULL;
 358         MatDuplicate(LS.getPetscMatrix(),MAT_COPY_VALUES,&_mat);
 359         _smb=NULL;
 360         VecDuplicate(LS.getPetscVector(),&_smb);                                ;
 361         _ksp=NULL;
 362         kspDuplicate(LS.getPetscKsp(),_mat,&_ksp);
 363         _prec=NULL;
 364         precDuplicate(LS.getPetscPc(),_ksp,&_prec);
 365         //    _ksp=LS.getPetscKsp();
 366         //    _prec=LS.getPetscPc();
 367         _isSparseMatrix=LS.isSparseMatrix();
 368 }
 369
 370 KSP
 371 LinearSolver::getPetscKsp() const
 372 {
 373         return (_ksp);
 374 }
 375
 376 Mat
 377 LinearSolver::getPetscMatrix() const
 378 {
 379         return (_mat);
 380 }
 381
 382 Vec
 383 LinearSolver::getPetscVector() const
 384 {
 385         return (_smb);
 386 }
 387
 388 void
 389 LinearSolver::viewPetscMatrix() const
 390 {
 391         MatView(_mat,PETSC_VIEWER_STDOUT_SELF);
 392 }
 393 void
 394 LinearSolver::viewPetscRHS() const
 395 {
 396         VecView(_smb,PETSC_VIEWER_STDOUT_SELF);
 397 }
 398 double
 399 LinearSolver::getPetscMatValue(int i, int j) const
 400 {
 401         double res;
 402         int idxm=i,idxn=j;
 403         MatGetValues(_mat,1,&idxm,1, &idxn,&res);
 404         return res;
 405 }
 406 double
 407 LinearSolver::getPetscRHSValue(int i) const
 408 {
 409         double res;
 410         int idxm=i;
 411         VecGetValues(_smb,1,&idxm,&res);
 412         return res;
 413 }
 414 PC
 415 LinearSolver::getPetscPc() const
 416 {
 417         return (_prec);
 418 }
 419
 420 Vector
 421 LinearSolver::solve( void )
 422 {
 423         if (_nameOfMethod.compare("GMRES")==0)
 424                 KSPSetType(_ksp,KSPGMRES);
 425         else if (_nameOfMethod.compare("LGMRES")==0)
 426                 KSPSetType(_ksp,KSPLGMRES);
 427         else if (_nameOfMethod.compare("CG")==0)
 428                 KSPSetType(_ksp,KSPCG);
 429         else if (_nameOfMethod.compare("BCGS")==0)
 430                 KSPSetType(_ksp,KSPBCGS);
 431         else if (_nameOfMethod.compare("CR")==0)
 432                 KSPSetType(_ksp,KSPCR);
 433         else if (_nameOfMethod.compare("CGS")==0)
 434                 KSPSetType(_ksp,KSPCGS);
 435         else if (_nameOfMethod.compare("BICG")==0)
 436                 KSPSetType(_ksp,KSPBICG);
 437         else if (_nameOfMethod.compare("GCR")==0)
 438                 KSPSetType(_ksp,KSPGCR);
 439         else if (_nameOfMethod.compare("LSQR")==0)
 440                 KSPSetType(_ksp,KSPLSQR);
 441         else if (_nameOfMethod.compare("CHOLESKY")==0)
 442         {
 443                 KSPSetType(_ksp,KSPGMRES);
 444                 PCSetType(_prec,PCCHOLESKY);
 445         }
 446         else if (_nameOfMethod.compare("LU")==0)
 447         {
 448                 KSPSetType(_ksp,KSPGMRES);
 449                 PCSetType(_prec,PCLU);
 450         }
 451         else
 452         {
 453                 string msg="Vector LinearSolver::solve( void ) : The method "+_nameOfMethod+" is not yet implemented.\n";
 454                 msg+="The methods implemented are : GMRES, BICG, CG, CHOLESKY, LU, BCGS, LGMRES, LSQR, CR, CGS and GCR.\n";
 455                 throw CdmathException(msg);
 456         }
 457
 458         if (_nameOfPc.compare("ILU")==0)
 459                 PCSetType(_prec,PCILU);
 460         else if (_nameOfPc.compare("LU")==0)
 461                 PCSetType(_prec,PCLU);
 462         else if (_nameOfPc.compare("ICC")==0)
 463                 PCSetType(_prec,PCICC);
 464         else if (_nameOfPc.compare("CHOLESKY")==0)
 465                 PCSetType(_prec,PCCHOLESKY);
 466         else if (_nameOfPc.compare("")==0)
 467                 PCSetType(_prec,PCNONE);
 468         else
 469         {
 470                 string msg="Vector LinearSolver::solve( void ) : The preconditioner "+_nameOfPc+" is not yet available.\n";
 471                 msg+="The preconditioners available are : void, ICC, ILU, CHOLESKY and LU.\n";
 472                 throw CdmathException(msg);
 473         }
 474
 475         KSPSetTolerances(_ksp,_tol*_tol*_tol*_tol,_tol,PETSC_DEFAULT,_numberMaxOfIter);
 476
 477         PetscInt its;
 478         PetscReal atol;
 479         PetscInt maxits;
 480
 481         Vec X;
 482         VecDuplicate(_smb,&X);
 483
 484         if (_isSingular)//Matrix should be symmetric semi-definite with constant vectors in its kernel
 485         {
 486                 //Check that the matrix is symmetric
 487                 PetscBool isSymetric;
 488                 MatIsSymmetric(_mat,_tol,&isSymetric);
 489                 if(!isSymetric)
 490                         {
 491                                 cout<<"Singular matrix is not symmetric, tolerance= "<< _tol<<endl;
 492                                 throw CdmathException("Singular matrix should be symmetric with kernel composed of constant vectors");
 493                         }
 494                 else
 495                 {
 496                         MatSetOption(_mat,MAT_HERMITIAN, PETSC_TRUE);
 497                         MatNullSpace nullsp;
 498                         MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);//Declaration of a kernel containing exclusively constant vectors
 499                         MatSetTransposeNullSpace(_mat, nullsp);//Transpose matrix has the same kernel since the matrix is symmetric
 500                         MatSetNullSpace(_mat, nullsp);
 501                         MatNullSpaceDestroy(&nullsp);
 502                 }
 503         }
 504
 505         if(_computeConditionNumber)
 506                 KSPSetComputeSingularValues(_ksp,PETSC_TRUE);
 507
 508         KSPSolve(_ksp,_smb,X);
 509
 510         KSPGetResidualNorm(_ksp,&atol);
 511         _residu=(double)atol;
 512
 513         KSPGetIterationNumber(_ksp,&its);
 514         _numberOfIter=(int)its;
 515
 516         KSPConvergedReason reason;
 517         KSPGetConvergedReason(_ksp,&reason);
 518
 519         if (reason==2 || reason==3)
 520                 _convergence=true;
 521         else{
 522                 _convergence=false;
 523                 cout<<"Linear system algorithm did not converge, divergence reason "<< reason <<endl;
 524         cout<<"Solver used "<<  _nameOfMethod<<", preconditioner "<<_nameOfPc<<endl;
 525                 cout<<"Final number of iteration= "<<_numberOfIter<<". Maximum allowed was " << _numberMaxOfIter<<endl;
 526                 cout<<"Final residual "<< _residu<< ". Objective was "<< _tol<<endl;
 527                 string msg="Linear system algorithm did not converge";
 528                 throw CdmathException(msg);
 529         }
 530
 531         if(_computeConditionNumber and !_isSingular)
 532         {
 533                 PetscReal sv_max, sv_min;
 534                 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
 535                 cout<<" Maximal singular value = " << sv_max <<", Minimal singular value = " << sv_min <<", Condition number = " << sv_max/sv_min <<endl;
 536         _conditionNumber=sv_max/sv_min;
 537         }
 538
 539         Vector X1=vecToVector(X);
 540
 541         return X1;
 542 }
 543
 544 Vec
 545 LinearSolver::vectorToVec(const Vector& myVector) const
 546 {
 547         int numberOfRows=myVector.getNumberOfRows();
 548         const double* values = myVector.getValues().getValues();
 549
 550         Vec result;
 551         VecCreate(PETSC_COMM_WORLD,&result);
 552         VecSetSizes(result,PETSC_DECIDE,numberOfRows);
 553         VecSetBlockSize(result,numberOfRows);
 554         VecSetFromOptions(result);
 555         int idx=0;//Index where to add the block of values
 556         VecSetValuesBlocked(result,1,&idx,values,INSERT_VALUES);
 557
 558         VecAssemblyBegin(result);
 559         VecAssemblyEnd(result);
 560
 561         return result;
 562 }
 563
 564 Vector
 565 LinearSolver::vecToVector(const Vec& vec) const
 566 {
 567         PetscInt numberOfRows;
 568         VecGetSize(vec,&numberOfRows);
 569     double * petscValues;
 570     VecGetArray(vec,&petscValues);
 571
 572     DoubleTab values (numberOfRows,petscValues);
 573         Vector result(numberOfRows);
 574     result.setValues(values);
 575
 576         return result;
 577 }
 578
 579 //----------------------------------------------------------------------
 580 const LinearSolver&
 581 LinearSolver::operator= ( const LinearSolver& linearSolver )
 582 //----------------------------------------------------------------------
 583 {
 584         _tol=linearSolver.getTolerance();
 585         _numberMaxOfIter=linearSolver.getNumberMaxOfIter();
 586         _nameOfMethod=linearSolver.getNameOfMethod();
 587         _residu=linearSolver.getResidu();
 588         _convergence=linearSolver.getStatus();
 589         _numberOfIter=linearSolver.getNumberOfIter();
 590         _isSingular=linearSolver.isMatrixSingular();
 591         _secondMember=linearSolver.getSndMember();
 592         _isSparseMatrix=linearSolver.isSparseMatrix();
 593         _nameOfPc="";
 594         setPreconditioner(linearSolver.getNameOfPc());
 595         _mat=NULL;
 596         MatDuplicate(linearSolver.getPetscMatrix(),MAT_COPY_VALUES,&_mat);
 597         _smb=NULL;
 598         VecDuplicate(linearSolver.getPetscVector(),&_smb);                              ;
 599         _ksp=NULL;
 600         kspDuplicate(linearSolver.getPetscKsp(),_mat,&_ksp);
 601         _prec=NULL;
 602         precDuplicate(linearSolver.getPetscPc(),_ksp,&_prec);
 603         return (*this);
 604 }