1 static char help[] = "Test MatMult() for Hermitian matrix.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc,char **args) 6 { 7 Mat A,As; 8 PetscBool flg; 9 PetscErrorCode ierr; 10 PetscMPIInt size; 11 PetscInt i,j; 12 PetscScalar v,sigma2; 13 PetscReal h2,sigma1=100.0; 14 PetscInt dim,Ii,J,n = 3,rstart,rend; 15 16 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 17 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 18 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL)); 19 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 20 dim = n*n; 21 22 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 23 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim)); 24 CHKERRQ(MatSetType(A,MATAIJ)); 25 CHKERRQ(MatSetFromOptions(A)); 26 CHKERRQ(MatSetUp(A)); 27 28 sigma2 = 10.0*PETSC_i; 29 h2 = 1.0/((n+1)*(n+1)); 30 31 CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 32 for (Ii=rstart; Ii<rend; Ii++) { 33 v = -1.0; i = Ii/n; j = Ii - i*n; 34 if (i>0) { 35 J = Ii-n; CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 36 } 37 if (i<n-1) { 38 J = Ii+n; CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 39 } 40 if (j>0) { 41 J = Ii-1; CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 42 } 43 if (j<n-1) { 44 J = Ii+1; CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 45 } 46 v = 4.0 - sigma1*h2; 47 CHKERRQ(MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES)); 48 } 49 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 50 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 51 52 /* Check whether A is symmetric */ 53 CHKERRQ(PetscOptionsHasName(NULL,NULL, "-check_symmetric", &flg)); 54 if (flg) { 55 CHKERRQ(MatIsSymmetric(A,0.0,&flg)); 56 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_USER,"A is not symmetric"); 57 } 58 CHKERRQ(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 59 60 /* make A complex Hermitian */ 61 Ii = 0; J = dim-1; 62 if (Ii >= rstart && Ii < rend) { 63 v = sigma2*h2; /* RealPart(v) = 0.0 */ 64 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 65 v = -sigma2*h2; 66 CHKERRQ(MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES)); 67 } 68 69 Ii = dim-2; J = dim-1; 70 if (Ii >= rstart && Ii < rend) { 71 v = sigma2*h2; /* RealPart(v) = 0.0 */ 72 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 73 v = -sigma2*h2; 74 CHKERRQ(MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES)); 75 } 76 77 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 78 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 79 CHKERRQ(MatViewFromOptions(A,NULL,"-disp_mat")); 80 81 /* Check whether A is Hermitian, then set A->hermitian flag */ 82 CHKERRQ(PetscOptionsHasName(NULL,NULL, "-check_Hermitian", &flg)); 83 if (flg && size == 1) { 84 CHKERRQ(MatIsHermitian(A,0.0,&flg)); 85 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_USER,"A is not Hermitian"); 86 } 87 CHKERRQ(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE)); 88 89 #if defined(PETSC_HAVE_SUPERLU_DIST) 90 /* Test Cholesky factorization */ 91 CHKERRQ(PetscOptionsHasName(NULL,NULL, "-test_choleskyfactor", &flg)); 92 if (flg) { 93 Mat F; 94 IS perm,iperm; 95 MatFactorInfo info; 96 PetscInt nneg,nzero,npos; 97 98 CHKERRQ(MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_CHOLESKY,&F)); 99 CHKERRQ(MatGetOrdering(A,MATORDERINGND,&perm,&iperm)); 100 CHKERRQ(MatCholeskyFactorSymbolic(F,A,perm,&info)); 101 CHKERRQ(MatCholeskyFactorNumeric(F,A,&info)); 102 103 CHKERRQ(MatGetInertia(F,&nneg,&nzero,&npos)); 104 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n",nneg,nzero,npos)); 105 CHKERRQ(MatDestroy(&F)); 106 CHKERRQ(ISDestroy(&perm)); 107 CHKERRQ(ISDestroy(&iperm)); 108 } 109 #endif 110 111 /* Create a Hermitian matrix As in sbaij format */ 112 CHKERRQ(MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&As)); 113 CHKERRQ(MatViewFromOptions(As,NULL,"-disp_mat")); 114 115 /* Test MatMult */ 116 CHKERRQ(MatMultEqual(A,As,10,&flg)); 117 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMult not equal"); 118 CHKERRQ(MatMultAddEqual(A,As,10,&flg)); 119 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatMultAdd not equal"); 120 121 /* Free spaces */ 122 CHKERRQ(MatDestroy(&A)); 123 CHKERRQ(MatDestroy(&As)); 124 ierr = PetscFinalize(); 125 return ierr; 126 } 127 128 /*TEST 129 130 build: 131 requires: complex 132 133 test: 134 args: -n 1000 135 output_file: output/ex127.out 136 137 test: 138 suffix: 2 139 nsize: 3 140 args: -n 1000 141 output_file: output/ex127.out 142 143 test: 144 suffix: superlu_dist 145 nsize: 3 146 requires: superlu_dist 147 args: -test_choleskyfactor -mat_superlu_dist_rowperm NOROWPERM 148 output_file: output/ex127_superlu_dist.out 149 TEST*/ 150