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