1 /* 2 Defines basic operations for the MATSEQAIJMKL matrix class. 3 This class is derived from the MATSEQAIJ class and retains the 4 compressed row storage (aka Yale sparse matrix format) but uses 5 sparse BLAS operations from the Intel Math Kernel Library (MKL) 6 wherever possible. 7 */ 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 #include <../src/mat/impls/aij/seq/aijmkl/aijmkl.h> 11 12 /* MKL include files. */ 13 #include <mkl_spblas.h> /* Sparse BLAS */ 14 15 typedef struct { 16 /* "Handle" used by SpMV2 inspector-executor routines. */ 17 sparse_matrix_t csrA; 18 } Mat_SeqAIJMKL; 19 20 extern PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat,MatAssemblyType); 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "MatConvert_SeqAIJMKL_SeqAIJ" 24 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJMKL_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) 25 { 26 /* This routine is only called to convert a MATAIJMKL to its base PETSc type, */ 27 /* so we will ignore 'MatType type'. */ 28 PetscErrorCode ierr; 29 Mat B = *newmat; 30 Mat_SeqAIJMKL *aijmkl=(Mat_SeqAIJMKL*)A->spptr; 31 32 PetscFunctionBegin; 33 if (reuse == MAT_INITIAL_MATRIX) { 34 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 35 } 36 37 /* Reset the original function pointers. */ 38 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ; 39 B->ops->destroy = MatDestroy_SeqAIJ; 40 B->ops->duplicate = MatDuplicate_SeqAIJ; 41 42 /* Free everything in the Mat_SeqAIJMKL data structure. 43 * We don't free the Mat_SeqAIJMKL struct itself, as this will 44 * cause problems later when MatDestroy() tries to free it. */ 45 /* Actually there is nothing to do here right now. 46 * When I've added use of the MKL SpMV2 inspector-executor routines, I should 47 * see if there is some way to clean up the "handle" used by SpMV2. */ 48 49 /* Change the type of B to MATSEQAIJ. */ 50 ierr = PetscObjectChangeTypeName((PetscObject)B, MATSEQAIJ);CHKERRQ(ierr); 51 52 *newmat = B; 53 PetscFunctionReturn(0); 54 } 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatDestroy_SeqAIJMKL" 58 PetscErrorCode MatDestroy_SeqAIJMKL(Mat A) 59 { 60 PetscErrorCode ierr; 61 Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL*) A->spptr; 62 63 PetscFunctionBegin; 64 /* Clean up everything in the Mat_SeqAIJMKL data structure, then free A->spptr. */ 65 mkl_sparse_destroy(aijmkl->csrA); 66 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 67 68 /* Change the type of A back to SEQAIJ and use MatDestroy_SeqAIJ() 69 * to destroy everything that remains. */ 70 ierr = PetscObjectChangeTypeName((PetscObject)A, MATSEQAIJ);CHKERRQ(ierr); 71 /* Note that I don't call MatSetType(). I believe this is because that 72 * is only to be called when *building* a matrix. I could be wrong, but 73 * that is how things work for the SuperLU matrix class. */ 74 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatDuplicate_SeqAIJMKL" 80 PetscErrorCode MatDuplicate_SeqAIJMKL(Mat A, MatDuplicateOption op, Mat *M) 81 { 82 PetscErrorCode ierr; 83 Mat_SeqAIJMKL *aijmkl = (Mat_SeqAIJMKL*) A->spptr; 84 Mat_SeqAIJMKL *aijmkl_dest = (Mat_SeqAIJMKL*) (*M)->spptr; 85 86 PetscFunctionBegin; 87 ierr = MatDuplicate_SeqAIJ(A,op,M);CHKERRQ(ierr); 88 ierr = PetscMemcpy((*M)->spptr,aijmkl,sizeof(Mat_SeqAIJMKL));CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatAssemblyEnd_SeqAIJMKL" 94 PetscErrorCode MatAssemblyEnd_SeqAIJMKL(Mat A, MatAssemblyType mode) 95 { 96 PetscErrorCode ierr; 97 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 98 99 PetscFunctionBegin; 100 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 101 102 /* Since a MATSEQAIJMKL matrix is really just a MATSEQAIJ with some 103 * extra information and some different methods, call the AssemblyEnd 104 * routine for a MATSEQAIJ. 105 * I'm not sure if this is the best way to do this, but it avoids 106 * a lot of code duplication. 107 * I also note that currently MATSEQAIJMKL doesn't know anything about 108 * the Mat_CompressedRow data structure that SeqAIJ now uses when there 109 * are many zero rows. If the SeqAIJ assembly end routine decides to use 110 * this, this may break things. (Don't know... haven't looked at it. 111 * Do I need to disable this somehow?) */ 112 a->inode.use = PETSC_FALSE; /* Must disable: otherwise the MKL routines won't get used. */ 113 ierr = MatAssemblyEnd_SeqAIJ(A, mode);CHKERRQ(ierr); 114 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "MatMult_SeqAIJMKL" 120 PetscErrorCode MatMult_SeqAIJMKL(Mat A,Vec xx,Vec yy) 121 { 122 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 123 const PetscScalar *x; 124 PetscScalar *y; 125 const MatScalar *aa; 126 PetscErrorCode ierr; 127 PetscInt m=A->rmap->n; 128 const PetscInt *aj,*ai; 129 PetscInt i; 130 131 /* Variables not in MatMult_SeqAIJ. */ 132 char transa = 'n'; /* Used to indicate to MKL that we are computing the transpose product. */ 133 134 PetscFunctionBegin; 135 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 136 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 137 aj = a->j; /* aj[k] gives column index for element aa[k]. */ 138 aa = a->a; /* Nonzero elements stored row-by-row. */ 139 ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */ 140 141 /* Call MKL sparse BLAS routine to do the MatMult. */ 142 mkl_cspblas_xcsrgemv(&transa,&m,aa,ai,aj,x,y); 143 144 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 145 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 146 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatMultAdd_SeqAIJMKL" 152 PetscErrorCode MatMultAdd_SeqAIJMKL(Mat A,Vec xx,Vec yy,Vec zz) 153 { 154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 155 const PetscScalar *x; 156 PetscScalar *y,*z; 157 const MatScalar *aa; 158 PetscErrorCode ierr; 159 PetscInt m=A->rmap->n; 160 const PetscInt *aj,*ai; 161 PetscInt i; 162 163 /* Variables not in MatMult_SeqAIJ. */ 164 char transa = 'n'; /* Used to indicate to MKL that we are computing the transpose product. */ 165 166 PetscFunctionBegin; 167 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 168 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 169 aj = a->j; /* aj[k] gives column index for element aa[k]. */ 170 aa = a->a; /* Nonzero elements stored row-by-row. */ 171 ai = a->i; /* ai[k] is the position in aa and aj where row k starts. */ 172 173 /* Call MKL sparse BLAS routine to do the MatMult. */ 174 mkl_cspblas_xcsrgemv(&transa,&m,aa,ai,aj,x,z); 175 176 /* Now add the vector y; MKL sparse BLAS does not have a MatMultAdd equivalent. */ 177 for (i=0; i<m; i++) { 178 z[i] += y[i]; 179 } 180 181 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 182 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 183 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 /* MatConvert_SeqAIJ_SeqAIJMKL converts a SeqAIJ matrix into a 188 * SeqAIJMKL matrix. This routine is called by the MatCreate_SeqMKLAIJ() 189 * routine, but can also be used to convert an assembled SeqAIJ matrix 190 * into a SeqAIJMKL one. */ 191 #undef __FUNCT__ 192 #define __FUNCT__ "MatConvert_SeqAIJ_SeqAIJMKL" 193 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJMKL(Mat A,MatType type,MatReuse reuse,Mat *newmat) 194 { 195 PetscErrorCode ierr; 196 Mat B = *newmat; 197 Mat_SeqAIJMKL *aijmkl; 198 199 PetscFunctionBegin; 200 if (reuse == MAT_INITIAL_MATRIX) { 201 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 202 } 203 204 ierr = PetscNewLog(B,&aijmkl);CHKERRQ(ierr); 205 B->spptr = (void*) aijmkl; 206 207 /* Set function pointers for methods that we inherit from AIJ but override. */ 208 B->ops->duplicate = MatDuplicate_SeqAIJMKL; 209 B->ops->assemblyend = MatAssemblyEnd_SeqAIJMKL; 210 B->ops->destroy = MatDestroy_SeqAIJMKL; 211 B->ops->mult = MatMult_SeqAIJMKL; 212 B->ops->multadd = MatMultAdd_SeqAIJMKL; 213 214 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaijmkl_seqaij_C",MatConvert_SeqAIJMKL_SeqAIJ);CHKERRQ(ierr); 215 216 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJMKL);CHKERRQ(ierr); 217 *newmat = B; 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "MatCreateSeqAIJMKL" 223 /*@C 224 MatCreateSeqAIJMKL - Creates a sparse matrix of type SEQAIJMKL. 225 This type inherits from AIJ and is largely identical, but uses sparse BLAS 226 routines from Intel MKL whenever possible. 227 Collective on MPI_Comm 228 229 Input Parameters: 230 + comm - MPI communicator, set to PETSC_COMM_SELF 231 . m - number of rows 232 . n - number of columns 233 . nz - number of nonzeros per row (same for all rows) 234 - nnz - array containing the number of nonzeros in the various rows 235 (possibly different for each row) or NULL 236 237 Output Parameter: 238 . A - the matrix 239 240 Notes: 241 If nnz is given then nz is ignored 242 243 Level: intermediate 244 245 .keywords: matrix, cray, sparse, parallel 246 247 .seealso: MatCreate(), MatCreateMPIAIJMKL(), MatSetValues() 248 @*/ 249 PetscErrorCode MatCreateSeqAIJMKL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 250 { 251 PetscErrorCode ierr; 252 253 PetscFunctionBegin; 254 ierr = MatCreate(comm,A);CHKERRQ(ierr); 255 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 256 ierr = MatSetType(*A,MATSEQAIJMKL);CHKERRQ(ierr); 257 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatCreate_SeqAIJMKL" 263 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJMKL(Mat A) 264 { 265 PetscErrorCode ierr; 266 267 PetscFunctionBegin; 268 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 269 ierr = MatConvert_SeqAIJ_SeqAIJMKL(A,MATSEQAIJMKL,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 274