CDMATH/linearsolver/src/SparseMatrixPetsc.cxx

   1 /*
   2  * SparseMatrixPetsc.cxx
   3  *
   4  *  Created on: 04/11/2017
   5  *      Author: mndjinga
   6  */
   7
   8 #include "SparseMatrixPetsc.hxx"
   9 #include "CdmathException.hxx"
  10
  11 using namespace std;
  12
  13 //----------------------------------------------------------------------
  14 SparseMatrixPetsc::SparseMatrixPetsc()
  15 //----------------------------------------------------------------------
  16 {
  17         _numberOfColumns=0;
  18         _numberOfRows=0;
  19         _numberOfNonZeros=0;
  20         _isSparseMatrix=true;
  21         _mat=NULL;
  22         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
  23 }
  24
  25 //----------------------------------------------------------------------
  26 SparseMatrixPetsc::SparseMatrixPetsc( int numberOfRows, int numberOfColumns)
  27 //----------------------------------------------------------------------
  28 {
  29         _numberOfRows = numberOfRows;
  30         _numberOfColumns=numberOfColumns;
  31         _isSparseMatrix=true;
  32         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
  33         MatCreateSeqAIJ(MPI_COMM_SELF,_numberOfRows,_numberOfColumns,PETSC_DEFAULT,NULL,&_mat);
  34 }
  35
  36 //----------------------------------------------------------------------
  37 SparseMatrixPetsc::SparseMatrixPetsc( Mat mat )
  38 //----------------------------------------------------------------------
  39 {
  40         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
  41         _isSparseMatrix=true;
  42         _mat=mat;
  43         //extract number of row and column
  44         MatGetSize(mat,&_numberOfRows,&_numberOfColumns);
  45
  46         //extract an upper bound for the total number of non zero coefficients
  47         MatInfo info;
  48         MatGetInfo(mat,MAT_LOCAL,&info);
  49         _numberOfNonZeros = info.nz_allocated;
  50 }
  51
  52 //----------------------------------------------------------------------
  53 SparseMatrixPetsc::SparseMatrixPetsc( int numberOfRows, int numberOfColumns, int nnz )
  54 //----------------------------------------------------------------------
  55 {
  56         _numberOfRows = numberOfRows;
  57         _numberOfColumns=numberOfColumns;
  58         _numberOfNonZeros=nnz;
  59         _isSparseMatrix=true;
  60         _mat=NULL;
  61         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
  62         MatCreateSeqAIJ(MPI_COMM_SELF,_numberOfRows,_numberOfColumns,_numberOfNonZeros,NULL,&_mat);
  63 }
  64
  65 //----------------------------------------------------------------------
  66 SparseMatrixPetsc::SparseMatrixPetsc( int blockSize, int numberOfRows, int numberOfColumns, int nnz )
  67 //----------------------------------------------------------------------
  68 {
  69         _numberOfRows = numberOfRows;
  70         _numberOfColumns=numberOfColumns;
  71         _numberOfNonZeros=nnz;
  72         _isSparseMatrix=true;
  73         _mat=NULL;
  74         PetscInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
  75         MatCreateSeqBAIJ(MPI_COMM_SELF,blockSize, _numberOfRows,_numberOfColumns,_numberOfNonZeros,NULL,&_mat);
  76 }
  77
  78 //----------------------------------------------------------------------
  79 SparseMatrixPetsc::SparseMatrixPetsc(const SparseMatrixPetsc& matrix)
  80 //----------------------------------------------------------------------
  81 {
  82         _isSparseMatrix=matrix.isSparseMatrix();
  83         MatDuplicate(matrix.getPetscMatrix(), MAT_COPY_VALUES,&_mat);
  84         MatGetSize(_mat,&_numberOfRows,&_numberOfColumns);
  85         //extract an upper bound for the total number of non zero coefficients
  86         MatInfo info;
  87         MatGetInfo(_mat,MAT_LOCAL,&info);
  88         _numberOfNonZeros = info.nz_allocated;
  89 }
  90
  91 SparseMatrixPetsc
  92 SparseMatrixPetsc::transpose() const
  93 {
  94         Mat mattranspose;
  95         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
  96         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
  97
  98         MatTranspose(_mat,MAT_INITIAL_MATRIX, &mattranspose);
  99         return SparseMatrixPetsc(mattranspose);
 100 }
 101
 102 void
 103 SparseMatrixPetsc::setValue( int i, int j, double value )
 104 {
 105         MatSetValues(_mat,1, &i, 1, &j, &value, INSERT_VALUES);
 106 }
 107
 108 void
 109 SparseMatrixPetsc::addValue( int i, int j, double value )
 110 {
 111         MatSetValues(_mat,1, &i, 1, &j, &value, ADD_VALUES);
 112 }
 113
 114 void
 115 SparseMatrixPetsc::setValue( int i, int j, Matrix M  )
 116 {
 117     int I,J;
 118     for (int k=0; k<M.getNumberOfRows(); k++)
 119         for (int l=0; l<M.getNumberOfColumns(); l++)
 120         {
 121             I=i+k;
 122             J=j+l;
 123             MatSetValues(_mat,1, &I, 1, &J, &M(k,l), INSERT_VALUES);
 124         }
 125 }
 126
 127 void
 128 SparseMatrixPetsc::addValue( int i, int j, Matrix M  )
 129 {
 130     int I,J;
 131     for (int k=0; k<M.getNumberOfRows(); k++)
 132         for (int l=0; l<M.getNumberOfColumns(); l++)
 133         {
 134             I=i+k;
 135             J=j+l;
 136             MatSetValues(_mat,1, &I, 1, &J, &M(k,l), ADD_VALUES);
 137         }
 138 }
 139
 140 void
 141 SparseMatrixPetsc::setValuesBlocked( int i, int j, Matrix M  )
 142 {
 143     int blockSize;
 144     MatGetBlockSize(_mat,&blockSize);
 145     if(blockSize!=M.getNumberOfRows() || blockSize!=M.getNumberOfColumns())
 146         throw CdmathException("SparseMatrixPetsc::setValuesBlocked : matrix size is different from sparse matrix block structure");
 147     double petscValues[blockSize*blockSize];
 148     for (int k=0; k<M.getNumberOfRows(); k++)
 149         for (int l=0; l<M.getNumberOfColumns(); l++)
 150             petscValues[k*blockSize+l]=M(k,l);
 151     MatSetValuesBlocked(_mat,1, &i, 1, &j, petscValues, INSERT_VALUES);
 152 }
 153
 154 void
 155 SparseMatrixPetsc::addValuesBlocked( int i, int j, Matrix M  )
 156 {
 157     int blockSize;
 158     MatGetBlockSize(_mat,&blockSize);
 159     if(blockSize!=M.getNumberOfRows() || blockSize!=M.getNumberOfColumns())
 160         throw CdmathException("SparseMatrixPetsc::addValuesBlocked : matrix size is different from sparse matrix block structure");
 161     double petscValues[blockSize*blockSize];
 162     for (int k=0; k<M.getNumberOfRows(); k++)
 163         for (int l=0; l<M.getNumberOfColumns(); l++)
 164             petscValues[k*blockSize+l]=M(k,l);
 165     MatSetValuesBlocked(_mat,1, &i, 1, &j, petscValues, ADD_VALUES);
 166 }
 167
 168 //----------------------------------------------------------------------
 169 double
 170 SparseMatrixPetsc::operator()( int i, int j ) const
 171 //----------------------------------------------------------------------
 172 {
 173         double res;
 174         int idxm=i,idxn=j;
 175         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 176         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 177
 178         MatGetValues(_mat,1,&idxm,1, &idxn,&res);
 179         return res;
 180 }
 181
 182 //----------------------------------------------------------------------
 183 SparseMatrixPetsc::~SparseMatrixPetsc()
 184 //----------------------------------------------------------------------
 185 {
 186         if(&_mat != NULL)
 187                 MatDestroy(&_mat);
 188         //PetscFinalize();
 189 }
 190
 191 Mat
 192 SparseMatrixPetsc::getPetscMatrix() const
 193 {
 194         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 195         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 196
 197         return (_mat);
 198 }
 199
 200 bool
 201 SparseMatrixPetsc::containsPetscMatrix() const
 202 {
 203         return true;
 204 }
 205 //----------------------------------------------------------------------
 206 const SparseMatrixPetsc&
 207 SparseMatrixPetsc::operator= ( const SparseMatrixPetsc& matrix )
 208 //----------------------------------------------------------------------
 209 {
 210         _isSparseMatrix=matrix.isSparseMatrix();
 211         MatDuplicate(matrix.getPetscMatrix(), MAT_COPY_VALUES,&_mat);
 212         MatGetSize(_mat,&_numberOfRows,&_numberOfColumns);
 213         //extract an upper bound for the total number of non zero coefficients
 214         MatInfo info;
 215         MatGetInfo(_mat,MAT_LOCAL,&info);
 216         _numberOfNonZeros = info.nz_allocated;
 217         return (*this);
 218 }
 219
 220 SparseMatrixPetsc
 221 operator+ (const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2)
 222 {
 223         int numberOfRows = matrix1.getNumberOfRows();
 224         int numberOfColumns = matrix1.getNumberOfColumns();
 225         int numberOfRows2 = matrix2.getNumberOfRows();
 226         int numberOfColumns2 = matrix2.getNumberOfColumns();
 227
 228         if(numberOfRows2!=numberOfRows || numberOfColumns2!=numberOfColumns)
 229         {
 230                 string msg="SparseMatrixPetsc::operator()+(const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2): number of rows or columns of the matrices is different!";
 231                 throw CdmathException(msg);
 232         }
 233
 234         Mat mat1=matrix1.getPetscMatrix();
 235         Mat mat2=matrix2.getPetscMatrix();
 236         Mat mat;
 237         MatDuplicate(mat1, MAT_COPY_VALUES,&mat);
 238         MatAXPY(mat,1,mat2,DIFFERENT_NONZERO_PATTERN);
 239
 240         return SparseMatrixPetsc(mat);
 241 }
 242
 243 SparseMatrixPetsc
 244 operator- (const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2)
 245 {
 246         int numberOfRows = matrix1.getNumberOfRows();
 247         int numberOfColumns = matrix1.getNumberOfColumns();
 248         int numberOfRows2 = matrix2.getNumberOfRows();
 249         int numberOfColumns2 = matrix2.getNumberOfColumns();
 250
 251         if(numberOfRows2!=numberOfRows || numberOfColumns2!=numberOfColumns)
 252         {
 253                 string msg="SparseMatrixPetsc::operator()-(const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2): number of rows or columns of the matrices is different!";
 254                 throw CdmathException(msg);
 255         }
 256         Mat mat1=matrix1.getPetscMatrix();
 257         Mat mat2=matrix2.getPetscMatrix();
 258         Mat mat;
 259         MatDuplicate(mat1, MAT_COPY_VALUES,&mat);
 260         MatAXPY(mat,-1,mat2,DIFFERENT_NONZERO_PATTERN);
 261
 262         return SparseMatrixPetsc(mat);
 263 }
 264
 265 SparseMatrixPetsc
 266 operator* (double value , const SparseMatrixPetsc& matrix )
 267 {
 268         Mat mat;
 269         MatDuplicate(matrix.getPetscMatrix(), MAT_COPY_VALUES,&mat);
 270         MatScale(mat, value);
 271
 272         return SparseMatrixPetsc(mat);
 273 }
 274
 275 SparseMatrixPetsc
 276 operator* (const SparseMatrixPetsc& matrix, double value )
 277 {
 278         Mat mat;
 279         MatDuplicate(matrix.getPetscMatrix(), MAT_COPY_VALUES,&mat);
 280         MatScale(mat, value);
 281
 282         return SparseMatrixPetsc(mat);
 283 }
 284
 285 SparseMatrixPetsc
 286 operator/ (const SparseMatrixPetsc& matrix, double value)
 287 {
 288         if(value==0.)
 289         {
 290                 string msg="SparseMatrixPetsc SparseMatrixPetsc::operator()/(const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2): division by zero";
 291                 throw CdmathException(msg);
 292         }
 293         Mat mat;
 294         MatDuplicate(matrix.getPetscMatrix(), MAT_COPY_VALUES,&mat);
 295         MatScale(mat, 1/value);
 296
 297         return SparseMatrixPetsc(mat);
 298 }
 299
 300 SparseMatrixPetsc
 301 operator*(const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2)
 302 {
 303         Mat mat1=matrix1.getPetscMatrix();
 304         Mat mat2=matrix2.getPetscMatrix();
 305         Mat mat;
 306         MatMatMult(mat1, mat2, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&mat);
 307
 308         return SparseMatrixPetsc(mat);
 309 }
 310
 311 SparseMatrixPetsc&
 312 SparseMatrixPetsc::operator*= (const SparseMatrixPetsc& matrix)
 313 {
 314         Mat mat1=matrix.getPetscMatrix();
 315         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 316         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 317
 318         Mat mat2;
 319         MatMatMult(_mat, mat1, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&mat2);
 320
 321         MatDestroy(&_mat);
 322         _mat=mat2;
 323         MatGetSize(_mat,&_numberOfRows,&_numberOfColumns);
 324         //extract an upper bound for the total number of non zero coefficients
 325         MatInfo info;
 326         MatGetInfo(_mat,MAT_LOCAL,&info);
 327         _numberOfNonZeros = info.nz_allocated;
 328         return (*this);
 329 }
 330
 331 Vector
 332 SparseMatrixPetsc::operator* (const Vector& vec) const
 333 {
 334         int numberOfRows=vec.getNumberOfRows();
 335         Vec X=vectorToVec(vec);
 336         Vec Y;
 337         VecDuplicate (X,&Y);
 338         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 339         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 340         MatMult(_mat,X,Y);
 341
 342         //Clean memory
 343         VecDestroy(&X);
 344
 345     Vector result=vecToVector(Y);
 346
 347         //Clean memory
 348         VecDestroy(&Y);
 349
 350         return result;
 351 }
 352
 353 SparseMatrixPetsc&
 354 SparseMatrixPetsc::operator+= (const SparseMatrixPetsc& matrix)
 355 {
 356         Mat mat1=matrix.getPetscMatrix();
 357         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 358         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 359         MatAXPY(_mat,1,mat1,DIFFERENT_NONZERO_PATTERN);
 360
 361         //extract an upper bound for the total number of non zero coefficients
 362         MatInfo info;
 363         MatGetInfo(_mat,MAT_LOCAL,&info);
 364         _numberOfNonZeros = info.nz_allocated;
 365
 366         return (*this);
 367 }
 368
 369 SparseMatrixPetsc&
 370 SparseMatrixPetsc::operator-= (const SparseMatrixPetsc& matrix)
 371 {
 372         Mat mat1=matrix.getPetscMatrix();
 373         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 374         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 375         MatAXPY(_mat,-1,mat1,DIFFERENT_NONZERO_PATTERN);
 376
 377         //extract an upper bound for the total number of non zero coefficients
 378         MatInfo info;
 379         MatGetInfo(_mat,MAT_LOCAL,&info);
 380         _numberOfNonZeros = info.nz_allocated;
 381
 382         return (*this);
 383 }
 384
 385 SparseMatrixPetsc&
 386 SparseMatrixPetsc::operator*= (double value)
 387 {
 388         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 389         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 390         MatScale(_mat, value);
 391         return (*this);
 392 }
 393
 394 SparseMatrixPetsc&
 395 SparseMatrixPetsc::operator/= (double value)
 396 {
 397         if(value==0.)
 398         {
 399                 string msg="SparseMatrixPetsc SparseMatrixPetsc::operator()/=(const SparseMatrixPetsc& matrix1, const SparseMatrixPetsc& matrix2): division by zero";
 400                 throw CdmathException(msg);
 401         }
 402         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 403         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 404         MatScale(_mat, 1/value);
 405         return (*this);
 406 }
 407
 408 void
 409 SparseMatrixPetsc::viewMatrix() const
 410 {
 411     MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 412         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 413
 414         MatView(_mat,PETSC_VIEWER_STDOUT_SELF);
 415 }
 416
 417 void
 418 SparseMatrixPetsc::saveMatrix(string filename, bool binaryMode) const
 419 {
 420     MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 421         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 422
 423         PetscViewer fileViewer;
 424         PetscViewerCreate(PETSC_COMM_WORLD,&fileViewer);
 425     PetscViewerFileSetMode(fileViewer,FILE_MODE_WRITE);
 426     PetscViewerFileSetName(fileViewer,filename.c_str());
 427
 428         if( binaryMode)
 429         {
 430                 PetscViewerSetType(fileViewer, PETSCVIEWERBINARY);
 431                 PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename.c_str(), &fileViewer);
 432         }
 433         else
 434         {
 435                 PetscViewerSetType(fileViewer, PETSCVIEWERASCII);
 436                 PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_WRITE, &fileViewer);
 437         }
 438
 439         MatView(_mat,fileViewer);
 440 }
 441
 442 double
 443 SparseMatrixPetsc::getMatrixCoeff(int i, int j) const
 444 {
 445         double res;
 446         int idxm=i,idxn=j;
 447         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 448         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 449
 450         MatGetValues(_mat,1,&idxm,1, &idxn,&res);
 451         return res;
 452 }
 453
 454 std::vector< double >
 455 SparseMatrixPetsc::getArray()
 456 {
 457         int size=_numberOfRows*_numberOfColumns;
 458
 459         vector< double >  result(size);
 460         double* values = result.data();
 461
 462         int * idxm = new int[_numberOfRows];
 463         int * idxn = new int[_numberOfColumns];
 464     for (int i=0;i<_numberOfRows;i++)
 465                 idxm[i]=i;
 466     for (int i=0;i<_numberOfColumns;i++)
 467                 idxn[i]=i;
 468
 469         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 470         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 471
 472         MatGetValues(_mat,_numberOfRows, idxm,_numberOfColumns, idxn,values);
 473
 474         return result;
 475 }
 476
 477 void
 478 SparseMatrixPetsc::diagonalShift(double lambda)
 479 {
 480     MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 481         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 482
 483     MatShift(_mat, lambda);
 484 }
 485
 486 void
 487 SparseMatrixPetsc::zeroEntries()
 488 {
 489     MatZeroEntries(_mat);
 490 }
 491
 492 Vector
 493 SparseMatrixPetsc::vecToVector(const Vec& vec) const
 494 {
 495         PetscInt numberOfRows;
 496         VecGetSize(vec,&numberOfRows);
 497     double * petscValues;
 498     VecGetArray(vec,&petscValues);
 499
 500     DoubleTab values (numberOfRows,petscValues);
 501         Vector result(numberOfRows);
 502     result.setValues(values);
 503
 504         return result;
 505 }
 506 Vec
 507 SparseMatrixPetsc::vectorToVec(const Vector& myVector) const
 508 {
 509         int numberOfRows=myVector.getNumberOfRows();
 510         const double* values = myVector.getValues().getValues();
 511
 512         Vec result;
 513         VecCreate(PETSC_COMM_WORLD,&result);
 514         VecSetSizes(result,PETSC_DECIDE,numberOfRows);
 515         VecSetBlockSize(result,numberOfRows);
 516         VecSetFromOptions(result);
 517         int idx=0;//Index where to add the block of values
 518         VecSetValuesBlocked(result,1,&idx,values,INSERT_VALUES);
 519
 520         VecAssemblyBegin(result);
 521         VecAssemblyEnd(result);
 522
 523         return result;
 524 }
 525
 526 bool SparseMatrixPetsc::isSymmetric(double tol) const
 527 {
 528         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
 529         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
 530
 531         //Check that the matrix is symmetric
 532         PetscBool isSymetric;
 533         MatIsSymmetric(_mat,tol,&isSymetric);
 534
 535         return isSymetric;
 536 }
 537
 538 int
 539 SparseMatrixPetsc::computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which, double tol) const
 540 {
 541   EPS            eps;         /* eigenproblem solver context */
 542   EPSType        type;
 543   PetscReal      error;
 544   PetscScalar    kr,ki;
 545   Vec            xr,xi;
 546   PetscInt       i,maxit,its, nconv;
 547
 548   SlepcInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
 549
 550   MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY);
 551   MatAssemblyEnd(_mat,MAT_FINAL_ASSEMBLY);
 552
 553   MatCreateVecs(_mat,&xr,NULL);
 554   MatCreateVecs(_mat,&xi,NULL);
 555
 556   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 557                 Create the eigensolver and set various options
 558      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 559   /*
 560      Create eigensolver context
 561   */
 562   EPSCreate(PETSC_COMM_WORLD,&eps);
 563
 564   /*
 565      Set operators. In this case, it is a standard eigenvalue problem
 566   */
 567   EPSSetOperators(eps,_mat,NULL);
 568   //if(isSymmetric(tol))
 569         //EPSSetProblemType(eps,EPS_HEP);
 570   //else
 571         EPSSetProblemType(eps,EPS_NHEP);
 572   EPSSetWhichEigenpairs(eps,which);
 573   EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
 574   EPSSetTolerances(eps,tol,PETSC_DEFAULT);
 575
 576   /*
 577      Set solver parameters at runtime
 578   */
 579   EPSSetFromOptions(eps);
 580
 581   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 582                       Solve the eigensystem
 583      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 584
 585   EPSSolve(eps);
 586   /*
 587      Optional: Get some information from the solver and display it
 588   */
 589   EPSGetIterationNumber(eps,&its);
 590   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
 591   EPSGetType(eps,&type);
 592   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 593   EPSGetDimensions(eps,&nev,NULL,NULL);
 594   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
 595   EPSGetTolerances(eps,&tol,&maxit);
 596   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
 597
 598   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 599                     Display solution and clean up
 600      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 601   /*
 602      Get number of converged approximate eigenpairs
 603   */
 604   EPSGetConverged(eps,&nconv);
 605   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
 606
 607   *valP=new double[nconv];
 608   *vecP=new double * [nconv];
 609   double * myVecp;
 610
 611   if (nconv>0) {
 612     /*
 613        Display eigenvalues and relative errors
 614     */
 615     PetscPrintf(PETSC_COMM_WORLD,
 616          "           k          ||Ax-kx||/||kx||\n"
 617          "   ----------------- ------------------\n");
 618
 619     for (int i=0;i<nconv;i++) {
 620       /*
 621         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
 622         ki (imaginary part)
 623       */
 624       EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
 625       /*
 626          Compute the relative error associated to each eigenpair
 627       */
 628       if(fabs(kr)>tol || fabs(ki)>tol)
 629       {
 630               EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);
 631
 632               if (ki!=0.0)
 633                 PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g\n",(double)kr,(double)ki,(double)error);
 634               else
 635                 PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",(double)kr,(double)error);
 636            }
 637            else
 638       {
 639               if (ki!=0.0)
 640                 PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12s\n",(double)kr,(double)ki,"Null eigenvalue");
 641               else
 642                 PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12s\n",(double)kr,"Null eigenvalue");
 643            }
 644
 645       *(*valP + i)=kr;
 646       VecGetArray(xr,&myVecp);
 647       *(*vecP+  i)=new double [_numberOfRows];
 648       memcpy(*(*vecP+  i),myVecp,_numberOfRows*sizeof(double)) ;
 649     }
 650     PetscPrintf(PETSC_COMM_WORLD,"\n");
 651     /*
 652      Free work space
 653     */
 654     EPSDestroy(&eps);
 655     VecDestroy(&xr);
 656     VecDestroy(&xi);
 657     SlepcFinalize();
 658
 659     return nconv;
 660   }
 661   else
 662   {
 663         /*
 664      Free work space
 665     */
 666     EPSDestroy(&eps);
 667     VecDestroy(&xr);
 668     VecDestroy(&xi);
 669     SlepcFinalize();
 670
 671     throw CdmathException("SparseMatrixPetsc::computeSpectrum : No eigenvector found");
 672   }
 673 }
 674
 675 int
 676 SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol) const
 677 {
 678   SVD            svd;         /* Singular value decomposition solver context */
 679   SVDType        type;
 680   PetscReal      error;
 681   PetscScalar    sigma;
 682   Vec            u,v;
 683   PetscInt       i,maxit,its, nconv;
 684
 685   SlepcInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
 686
 687   MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY);
 688   MatAssemblyEnd(  _mat,MAT_FINAL_ASSEMBLY);
 689
 690   MatCreateVecs(_mat,&u,NULL);
 691   MatCreateVecs(_mat,NULL,&v);
 692
 693   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 694                 Create the singular value solver and set various options
 695      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 696   /*
 697      Create singular value decomposition context
 698   */
 699   SVDCreate(PETSC_COMM_WORLD,&svd);
 700
 701   /*
 702      Set operators. In this case, it is a standard singular value problem
 703   */
 704 #if (SLEPC_VERSION_MAJOR==3) && (SLEPC_VERSION_MINOR > 14)
 705         SVDSetOperators(svd,_mat,NULL);
 706 #else
 707         SVDSetOperator(svd,_mat);
 708 #endif
 709
 710   SVDSetWhichSingularTriplets(svd,which);
 711   SVDSetDimensions(svd,nsv,PETSC_DEFAULT,PETSC_DEFAULT);
 712   SVDSetTolerances(svd,tol,PETSC_DEFAULT);
 713
 714   /*
 715      Set solver parameters at runtime
 716   */
 717   SVDSetFromOptions(svd);
 718
 719   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 720                       Solve the SVD problem
 721      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 722
 723   SVDSolve(svd);
 724   /*
 725      Optional: Get some information from the solver and display it
 726   */
 727   SVDGetIterationNumber(svd,&its);
 728   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
 729   SVDGetType(svd,&type);
 730   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 731   SVDGetDimensions(svd,&nsv,NULL,NULL);
 732   PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %D\n",nsv);
 733   SVDGetTolerances(svd,&tol,&maxit);
 734   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
 735
 736   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 737                     Display solution and clean up
 738      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 739   /*
 740      Get number of converged approximate singular values
 741   */
 742   SVDGetConverged(svd,&nconv);
 743   PetscPrintf(PETSC_COMM_WORLD," Number of converged singular values: %D\n\n",nconv);
 744
 745   *valS=new double[nconv];
 746   *vecS=new double * [2*nconv];
 747   double * myVecS;
 748
 749   if (nconv>0) {
 750     /*
 751        Display eigenvalues and relative errors
 752     */
 753     PetscPrintf(PETSC_COMM_WORLD,
 754          "           k          ||Ax-kx||/||kx||\n"
 755          "   ----------------- ------------------\n");
 756
 757     for (int i=0;i<nconv;i++) {
 758       /*
 759         Get converged singular values: i-th eigenvalue is stored in valS
 760       */
 761       SVDGetSingularTriplet(svd,i,*valS+i, u, v);
 762       if(fabs(*(*valS+i))>tol)
 763       {
 764               /*
 765                  Compute the relative error associated to each singular value
 766               */
 767               SVDComputeError( svd, i, SVD_ERROR_RELATIVE, &error );
 768
 769               PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",(double)*(*valS+i),(double)error);
 770                 }
 771        else
 772           PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12s\n",(double)*(*valS+i),"Null singular value");
 773
 774       VecGetArray(u,&myVecS);
 775       *(*vecS + i)=new double [_numberOfRows];
 776       memcpy(*(*vecS + i),myVecS,_numberOfRows*sizeof(double)) ;
 777       VecGetArray(v,&myVecS);
 778       *(*vecS + i + nconv)=new double [_numberOfColumns];
 779       memcpy(*(*vecS + i + nconv),myVecS,_numberOfColumns*sizeof(double)) ;
 780     }
 781     PetscPrintf(PETSC_COMM_WORLD,"\n");
 782     /*
 783      Free work space
 784     */
 785     SVDDestroy(&svd);
 786     VecDestroy(&u);
 787     VecDestroy(&v);
 788     SlepcFinalize();
 789
 790     return nconv;
 791   }
 792   else
 793   {
 794         /*
 795      Free work space
 796     */
 797     SVDDestroy(&svd);
 798     VecDestroy(&u);
 799     VecDestroy(&v);
 800     SlepcFinalize();
 801
 802     throw CdmathException("SparseMatrixPetsc::computeSVD : singular value decomposition failed");
 803   }
 804 }
 805
 806 std::vector< double >
 807 SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol) const
 808 {
 809         int nconv;
 810         double * valP;
 811         double **vecP;
 812
 813         nconv=computeSpectrum(nev, &valP, &vecP, which, tol);
 814
 815     std::vector< double > result(nconv);
 816
 817     for (int i=0;i<nconv;i++)
 818         result[i]=valP[i];
 819
 820         delete[] valP;
 821     for (int i=0;i<nconv;i++)
 822                 delete[] vecP[i];
 823         delete[] vecP;
 824
 825     return result;
 826 }
 827
 828 std::vector< Vector >
 829 SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol) const
 830 {
 831         int nconv;
 832         double * valP;
 833         double **vecP;
 834
 835         nconv=computeSpectrum(nev, &valP, &vecP, which, tol);
 836
 837     std::vector< Vector > result(nconv);
 838
 839     for (int i=0;i<nconv;i++)
 840     {
 841                 DoubleTab values (_numberOfRows,vecP[i]);
 842         Vector myVecP(_numberOfRows);
 843         myVecP.setValues(values);
 844         result[i]=myVecP;
 845         }
 846
 847         delete[] valP;
 848     for (int i=0;i<nconv;i++)
 849                 delete[] vecP[i];
 850         delete[] vecP;
 851
 852     return result;
 853 }
 854
 855 MEDCoupling::DataArrayDouble *
 856 SparseMatrixPetsc::getEigenvectorsDataArrayDouble(int nev, EPSWhich which, double tol) const
 857 {
 858         int nconv;
 859         double * valP;
 860         double **vecP;
 861
 862         nconv=computeSpectrum(nev, &valP, &vecP, which, tol);
 863
 864 #ifdef MEDCoupling_VERSION_VERSION_GREATER_9_4
 865         std::vector< long unsigned int > compoId(1);
 866 #else
 867         std::vector< int > compoId(1);
 868 #endif
 869         MEDCoupling::DataArrayDouble *arrays=MEDCoupling::DataArrayDouble::New();
 870         MEDCoupling::DataArrayDouble *array =MEDCoupling::DataArrayDouble::New();
 871         arrays->alloc(_numberOfRows,nconv);
 872
 873     for (int i=0;i<nconv;i++)
 874     {
 875                 array->useArray(vecP[i],true, MEDCoupling::DeallocType::CPP_DEALLOC, _numberOfRows,1);
 876                 compoId[0]=i;
 877                 arrays->setSelectedComponents(array,compoId);
 878                 arrays->setInfoOnComponent(i,std::to_string(valP[i]));
 879         }
 880         delete[] valP;
 881         delete[] vecP;
 882
 883     return arrays;
 884 }
 885
 886 double
 887 SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const
 888 {
 889         if(isSymmetric(tol))//Eigenvalues computation is more robust that singular values computation
 890         {
 891                 int nev, nconv;
 892                 double lambda_max, lambda_min;
 893                 std::vector< double > my_ev;
 894
 895                 /*** Lowest eigenvalue, first check if matrix is singular ***/
 896             if(isSingular)
 897                 nev=2;
 898             else
 899                         nev=1 ;
 900
 901                 my_ev=getEigenvalues( nev, EPS_SMALLEST_MAGNITUDE, tol);
 902                 nconv=my_ev.size();
 903                 if(nconv<nev)
 904                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the smallest eigenvalue");
 905             lambda_min=my_ev[nev-1];
 906
 907             /*** Largest eigenvalue ***/
 908             nev=1;
 909                 my_ev=getEigenvalues( nev, EPS_LARGEST_MAGNITUDE, tol);
 910                 nconv=my_ev.size();
 911                 if(nconv<nev)
 912                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the largest eigenvalue");
 913             lambda_max=my_ev[nev-1];
 914
 915             return fabs(lambda_max)/fabs(lambda_min);
 916         }
 917         else//Singular values computation is more robust than eigenvalues computation
 918         {
 919                 int nsv, nconv;
 920                 double * valS;
 921                 double ** vecS;
 922                 double sigma_max, sigma_min;
 923
 924                 /*** Lowest singular value, first check if matrix is singular ***/
 925             if(isSingular)
 926                 nsv=2;
 927             else
 928                         nsv=1 ;
 929
 930                 nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol);
 931                 if(nconv<nsv)
 932                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the smallest singular value");
 933             sigma_min=valS[nsv-1];
 934             delete[] valS, vecS;
 935
 936             /*** Largest singular value ***/
 937             nsv=1;
 938                 nconv=computeSVD(nsv, &valS, &vecS, SVD_LARGEST, tol);
 939                 if(nconv<nsv)
 940                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the largest singular value");
 941             sigma_max=valS[nsv-1];
 942             delete[] valS, vecS;
 943
 944             return sigma_max/sigma_min;
 945         }
 946 }
 947
 948 std::vector< double >
 949 SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol) const
 950 {
 951         int nconv;
 952         double * valS;
 953         double **vecS;
 954
 955         nconv=computeSVD(nsv, &valS, &vecS, which, tol);
 956
 957     std::vector< double > result(nconv);
 958
 959     for (int i=0;i<nconv;i++)
 960         result[i]=valS[i];
 961
 962         delete[] valS, vecS;
 963
 964     return result;
 965 }
 966
 967 std::vector< Vector >
 968 SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol) const
 969 {
 970         int nconv;
 971         double * valS;
 972         double **vecS;
 973
 974         nconv=computeSVD(nsv, &valS, &vecS, which, tol);
 975
 976     std::vector< Vector > result(2*nconv);
 977
 978     for (int i=0;i<nconv;i++)
 979     {
 980                 DoubleTab values_left (_numberOfRows,vecS[i]);
 981         Vector myVecS_left(_numberOfRows);
 982         myVecS_left.setValues(values_left);
 983         result[i]=myVecS_left;
 984                 DoubleTab values_right (_numberOfColumns,vecS[i+nconv]);
 985         Vector myVecS_right(_numberOfColumns);
 986         myVecS_right.setValues(values_right);
 987         result[i+nconv]=myVecS_right;
 988         }
 989
 990         delete[] valS;
 991     for (int i=0;i<nconv;i++)
 992                 delete[] vecS[i];
 993         delete[] vecS;
 994
 995     return result;
 996 }
 997
 998 void
 999 SparseMatrixPetsc::leftDiagonalScale(Vector v)
1000 {
1001         Vec vec=vectorToVec(v);
1002
1003         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
1004         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
1005
1006         MatDiagonalScale(_mat, vec, NULL);
1007 }
1008 void
1009 SparseMatrixPetsc::rightDiagonalScale(Vector v)
1010 {
1011         Vec vec=vectorToVec(v);
1012
1013         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
1014         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
1015
1016         MatDiagonalScale(_mat, NULL, vec);
1017 }