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