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, EPSType type) const
 540 {
 541   EPS            eps;         /* eigenproblem solver context */
 542   PetscReal      error;
 543   PetscScalar    kr,ki;
 544   Vec            xr,xi;
 545   PetscInt       i,maxit,its, nconv;
 546
 547   SlepcInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
 548
 549   MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY);
 550   MatAssemblyEnd(_mat,MAT_FINAL_ASSEMBLY);
 551
 552   MatCreateVecs(_mat,&xr,NULL);
 553   MatCreateVecs(_mat,&xi,NULL);
 554
 555   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 556                 Create the eigensolver and set various options
 557      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 558   /*
 559      Create eigensolver context
 560   */
 561   EPSCreate(PETSC_COMM_WORLD,&eps);
 562
 563   /*
 564      Set operators. In this case, it is a standard eigenvalue problem
 565   */
 566   EPSSetOperators(eps,_mat,NULL);
 567   //if(isSymmetric(tol))
 568         //EPSSetProblemType(eps,EPS_HEP);
 569   //else
 570         EPSSetProblemType(eps,EPS_NHEP);
 571   EPSSetWhichEigenpairs(eps,which);
 572   EPSSetDimensions(eps,nev,PETSC_DEFAULT,PETSC_DEFAULT);
 573   EPSSetTolerances(eps,tol,PETSC_DEFAULT);
 574   EPSSetType(eps,type);
 575   /*
 576      Set solver parameters at runtime
 577   */
 578   EPSSetFromOptions(eps);
 579
 580   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 581                       Solve the eigensystem
 582      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 583
 584   EPSSolve(eps);
 585   /*
 586      Optional: Get some information from the solver and display it
 587   */
 588   EPSGetIterationNumber(eps,&its);
 589   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
 590   EPSGetType(eps,&type);
 591   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 592   EPSGetDimensions(eps,&nev,NULL,NULL);
 593   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
 594   EPSGetTolerances(eps,&tol,&maxit);
 595   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
 596
 597   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 598                     Display solution and clean up
 599      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 600   /*
 601      Get number of converged approximate eigenpairs
 602   */
 603   EPSGetConverged(eps,&nconv);
 604   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
 605
 606   *valP=new double[nconv];
 607   *vecP=new double * [nconv];
 608   double * myVecp;
 609
 610   if (nconv>0) {
 611     /*
 612        Display eigenvalues and relative errors
 613     */
 614     PetscPrintf(PETSC_COMM_WORLD,
 615          "           k          ||Ax-kx||/||kx||\n"
 616          "   ----------------- ------------------\n");
 617
 618     for (int i=0;i<nconv;i++) {
 619       /*
 620         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
 621         ki (imaginary part)
 622       */
 623       EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
 624       /*
 625          Compute the relative error associated to each eigenpair
 626       */
 627       if(fabs(kr)>tol || fabs(ki)>tol)
 628       {
 629               EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);
 630
 631               if (ki!=0.0)
 632                 PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g\n",(double)kr,(double)ki,(double)error);
 633               else
 634                 PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",(double)kr,(double)error);
 635            }
 636            else
 637       {
 638               if (ki!=0.0)
 639                 PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12s\n",(double)kr,(double)ki,"Null eigenvalue");
 640               else
 641                 PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12s\n",(double)kr,"Null eigenvalue");
 642            }
 643
 644       *(*valP + i)=kr;
 645       VecGetArray(xr,&myVecp);
 646       *(*vecP+  i)=new double [_numberOfRows];
 647       memcpy(*(*vecP+  i),myVecp,_numberOfRows*sizeof(double)) ;
 648     }
 649     PetscPrintf(PETSC_COMM_WORLD,"\n");
 650     /*
 651      Free work space
 652     */
 653     EPSDestroy(&eps);
 654     VecDestroy(&xr);
 655     VecDestroy(&xi);
 656     SlepcFinalize();
 657
 658     return nconv;
 659   }
 660   else
 661   {
 662         /*
 663      Free work space
 664     */
 665     EPSDestroy(&eps);
 666     VecDestroy(&xr);
 667     VecDestroy(&xi);
 668     SlepcFinalize();
 669
 670     throw CdmathException("SparseMatrixPetsc::computeSpectrum : No eigenvector found");
 671   }
 672 }
 673
 674 int
 675 SparseMatrixPetsc::computeSVD(int nsv, double ** valS, double ***vecS, SVDWhich which, double tol, SVDType type) const
 676 {
 677   SVD            svd;         /* Singular value decomposition solver context */
 678   PetscReal      error;
 679   PetscScalar    sigma;
 680   Vec            u,v;
 681   PetscInt       i,maxit,its, nconv;
 682
 683   SlepcInitialize(0, (char ***)"", PETSC_NULL, PETSC_NULL);
 684
 685   MatAssemblyBegin(_mat,MAT_FINAL_ASSEMBLY);
 686   MatAssemblyEnd(  _mat,MAT_FINAL_ASSEMBLY);
 687
 688   MatCreateVecs(_mat,&u,NULL);
 689   MatCreateVecs(_mat,NULL,&v);
 690
 691   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 692                 Create the singular value solver and set various options
 693      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 694   /*
 695      Create singular value decomposition context
 696   */
 697   SVDCreate(PETSC_COMM_WORLD,&svd);
 698
 699   /*
 700      Set operators. In this case, it is a standard singular value problem
 701   */
 702 #if (SLEPC_VERSION_MAJOR==3) && (SLEPC_VERSION_MINOR > 14)
 703         SVDSetOperators(svd,_mat,NULL);
 704 #else
 705         SVDSetOperator(svd,_mat);
 706 #endif
 707
 708   SVDSetWhichSingularTriplets(svd,which);
 709   SVDSetDimensions(svd,nsv,PETSC_DEFAULT,PETSC_DEFAULT);
 710   SVDSetType(svd, type);
 711   SVDSetTolerances(svd,tol,PETSC_DEFAULT);
 712
 713   /*
 714      Set solver parameters at runtime
 715   */
 716   SVDSetFromOptions(svd);
 717
 718   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 719                       Solve the SVD problem
 720      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 721
 722   SVDSolve(svd);
 723   /*
 724      Optional: Get some information from the solver and display it
 725   */
 726   SVDGetIterationNumber(svd,&its);
 727   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
 728   SVDGetType(svd,&type);
 729   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 730   SVDGetDimensions(svd,&nsv,NULL,NULL);
 731   PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %D\n",nsv);
 732   SVDGetTolerances(svd,&tol,&maxit);
 733   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
 734
 735   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 736                     Display solution and clean up
 737      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 738   /*
 739      Get number of converged approximate singular values
 740   */
 741   SVDGetConverged(svd,&nconv);
 742   PetscPrintf(PETSC_COMM_WORLD," Number of converged singular values: %D\n\n",nconv);
 743
 744   *valS=new double[nconv];
 745   *vecS=new double * [2*nconv];
 746   double * myVecS;
 747
 748   if (nconv>0) {
 749     /*
 750        Display eigenvalues and relative errors
 751     */
 752     PetscPrintf(PETSC_COMM_WORLD,
 753          "           k          ||Ax-kx||/||kx||\n"
 754          "   ----------------- ------------------\n");
 755
 756     for (int i=0;i<nconv;i++) {
 757       /*
 758         Get converged singular values: i-th eigenvalue is stored in valS
 759       */
 760       SVDGetSingularTriplet(svd,i,*valS+i, u, v);
 761       if(fabs(*(*valS+i))>tol)
 762       {
 763               /*
 764                  Compute the relative error associated to each singular value
 765               */
 766               SVDComputeError( svd, i, SVD_ERROR_RELATIVE, &error );
 767
 768               PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",(double)*(*valS+i),(double)error);
 769                 }
 770        else
 771           PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12s\n",(double)*(*valS+i),"Null singular value");
 772
 773       VecGetArray(u,&myVecS);
 774       *(*vecS + i)=new double [_numberOfRows];
 775       memcpy(*(*vecS + i),myVecS,_numberOfRows*sizeof(double)) ;
 776       VecGetArray(v,&myVecS);
 777       *(*vecS + i + nconv)=new double [_numberOfColumns];
 778       memcpy(*(*vecS + i + nconv),myVecS,_numberOfColumns*sizeof(double)) ;
 779     }
 780     PetscPrintf(PETSC_COMM_WORLD,"\n");
 781     /*
 782      Free work space
 783     */
 784     SVDDestroy(&svd);
 785     VecDestroy(&u);
 786     VecDestroy(&v);
 787     SlepcFinalize();
 788
 789     return nconv;
 790   }
 791   else
 792   {
 793         /*
 794      Free work space
 795     */
 796     SVDDestroy(&svd);
 797     VecDestroy(&u);
 798     VecDestroy(&v);
 799     SlepcFinalize();
 800
 801     throw CdmathException("SparseMatrixPetsc::computeSVD : singular value decomposition failed");
 802   }
 803 }
 804
 805 std::vector< double >
 806 SparseMatrixPetsc::getEigenvalues(int nev, EPSWhich which, double tol, EPSType type) const
 807 {
 808         int nconv;
 809         double * valP;
 810         double **vecP;
 811
 812         nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type);
 813
 814     std::vector< double > result(nconv);
 815
 816     for (int i=0;i<nconv;i++)
 817         result[i]=valP[i];
 818
 819         delete[] valP;
 820     for (int i=0;i<nconv;i++)
 821                 delete[] vecP[i];
 822         delete[] vecP;
 823
 824     return result;
 825 }
 826
 827 std::vector< Vector >
 828 SparseMatrixPetsc::getEigenvectors(int nev, EPSWhich which, double tol, EPSType type) const
 829 {
 830         int nconv;
 831         double * valP;
 832         double **vecP;
 833
 834         nconv=computeSpectrum(nev, &valP, &vecP, which, tol,type);
 835
 836     std::vector< Vector > result(nconv);
 837
 838     for (int i=0;i<nconv;i++)
 839     {
 840                 DoubleTab values (_numberOfRows,vecP[i]);
 841         Vector myVecP(_numberOfRows);
 842         myVecP.setValues(values);
 843         result[i]=myVecP;
 844         }
 845
 846         delete[] valP;
 847     for (int i=0;i<nconv;i++)
 848                 delete[] vecP[i];
 849         delete[] vecP;
 850
 851     return result;
 852 }
 853
 854 MEDCoupling::DataArrayDouble *
 855 SparseMatrixPetsc::getEigenvectorsDataArrayDouble(int nev, EPSWhich which, double tol, EPSType type) const
 856 {
 857         int nconv;
 858         double * valP;
 859         double **vecP;
 860
 861         nconv=computeSpectrum(nev, &valP, &vecP, which, tol);
 862
 863 #ifdef MEDCoupling_VERSION_VERSION_GREATER_9_4
 864         std::vector< long unsigned int > compoId(1);
 865 #else
 866         std::vector< int > compoId(1);
 867 #endif
 868         MEDCoupling::DataArrayDouble *arrays=MEDCoupling::DataArrayDouble::New();
 869         MEDCoupling::DataArrayDouble *array =MEDCoupling::DataArrayDouble::New();
 870         arrays->alloc(_numberOfRows,nconv);
 871
 872     for (int i=0;i<nconv;i++)
 873     {
 874                 array->useArray(vecP[i],true, MEDCoupling::DeallocType::CPP_DEALLOC, _numberOfRows,1);
 875                 compoId[0]=i;
 876                 arrays->setSelectedComponents(array,compoId);
 877                 arrays->setInfoOnComponent(i,std::to_string(valP[i]));
 878         }
 879         delete[] valP;
 880         delete[] vecP;
 881
 882     return arrays;
 883 }
 884
 885 double
 886 SparseMatrixPetsc::getConditionNumber(bool isSingular, double tol) const
 887 {
 888         if(isSymmetric(tol))//Eigenvalues computation is more robust that singular values computation
 889         {
 890                 int nev, nconv;
 891                 double lambda_max, lambda_min;
 892                 std::vector< double > my_ev;
 893
 894                 /*** Lowest eigenvalue, first check if matrix is singular ***/
 895             if(isSingular)
 896                 nev=2;
 897             else
 898                         nev=1 ;
 899
 900                 my_ev=getEigenvalues( nev, EPS_SMALLEST_MAGNITUDE, tol);
 901                 nconv=my_ev.size();
 902                 if(nconv<nev)
 903                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the smallest eigenvalue");
 904             lambda_min=my_ev[nev-1];
 905
 906             /*** Largest eigenvalue ***/
 907             nev=1;
 908                 my_ev=getEigenvalues( nev, EPS_LARGEST_MAGNITUDE, tol);
 909                 nconv=my_ev.size();
 910                 if(nconv<nev)
 911                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the largest eigenvalue");
 912             lambda_max=my_ev[nev-1];
 913
 914             return fabs(lambda_max)/fabs(lambda_min);
 915         }
 916         else//Singular values computation is more robust than eigenvalues computation
 917         {
 918                 int nsv, nconv;
 919                 double * valS;
 920                 double ** vecS;
 921                 double sigma_max, sigma_min;
 922
 923                 /*** Lowest singular value, first check if matrix is singular ***/
 924             if(isSingular)
 925                 nsv=2;
 926             else
 927                         nsv=1 ;
 928
 929                 nconv=computeSVD(nsv, &valS, &vecS, SVD_SMALLEST, tol,SVDCYCLIC);
 930                 if(nconv<nsv)
 931                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the smallest singular value");
 932             sigma_min=valS[nsv-1];
 933             delete[] valS, vecS;
 934
 935             /*** Largest singular value ***/
 936             nsv=1;
 937                 nconv=computeSVD(nsv, &valS, &vecS, SVD_LARGEST, tol,SVDCYCLIC);
 938                 if(nconv<nsv)
 939                         throw CdmathException("SparseMatrixPetsc::getConditionNumber could not find the largest singular value");
 940             sigma_max=valS[nsv-1];
 941             delete[] valS, vecS;
 942
 943             return sigma_max/sigma_min;
 944         }
 945 }
 946
 947 std::vector< double >
 948 SparseMatrixPetsc::getSingularValues(int nsv, SVDWhich which, double tol, SVDType type) const
 949 {
 950         int nconv;
 951         double * valS;
 952         double **vecS;
 953
 954         nconv=computeSVD(nsv, &valS, &vecS, which, tol, type);
 955
 956     std::vector< double > result(nconv);
 957
 958     for (int i=0;i<nconv;i++)
 959         result[i]=valS[i];
 960
 961         delete[] valS, vecS;
 962
 963     return result;
 964 }
 965
 966 std::vector< Vector >
 967 SparseMatrixPetsc::getSingularVectors(int nsv, SVDWhich which, double tol, SVDType type) const
 968 {
 969         int nconv;
 970         double * valS;
 971         double **vecS;
 972
 973         nconv=computeSVD(nsv, &valS, &vecS, which, tol, type);
 974
 975     std::vector< Vector > result(2*nconv);
 976
 977     for (int i=0;i<nconv;i++)
 978     {
 979                 DoubleTab values_left (_numberOfRows,vecS[i]);
 980         Vector myVecS_left(_numberOfRows);
 981         myVecS_left.setValues(values_left);
 982         result[i]=myVecS_left;
 983                 DoubleTab values_right (_numberOfColumns,vecS[i+nconv]);
 984         Vector myVecS_right(_numberOfColumns);
 985         myVecS_right.setValues(values_right);
 986         result[i+nconv]=myVecS_right;
 987         }
 988
 989         delete[] valS;
 990     for (int i=0;i<nconv;i++)
 991                 delete[] vecS[i];
 992         delete[] vecS;
 993
 994     return result;
 995 }
 996
 997 void
 998 SparseMatrixPetsc::leftDiagonalScale(Vector v)
 999 {
1000         Vec vec=vectorToVec(v);
1001
1002         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
1003         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
1004
1005         MatDiagonalScale(_mat, vec, NULL);
1006 }
1007 void
1008 SparseMatrixPetsc::rightDiagonalScale(Vector v)
1009 {
1010         Vec vec=vectorToVec(v);
1011
1012         MatAssemblyBegin(_mat, MAT_FINAL_ASSEMBLY);
1013         MatAssemblyEnd(_mat, MAT_FINAL_ASSEMBLY);
1014
1015         MatDiagonalScale(_mat, NULL, vec);
1016 }