1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatConvert_SeqDense_SeqAIJ" 13 PETSC_EXTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 14 { 15 Mat B; 16 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 17 PetscErrorCode ierr; 18 PetscInt i, j; 19 PetscInt *rows, *nnz; 20 MatScalar *aa = a->v, *vals; 21 22 PetscFunctionBegin; 23 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 24 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 25 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 26 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 27 for (j=0; j<A->cmap->n; j++) { 28 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i]; 29 aa += a->lda; 30 } 31 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 32 aa = a->v; 33 for (j=0; j<A->cmap->n; j++) { 34 PetscInt numRows = 0; 35 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];} 36 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 37 aa += a->lda; 38 } 39 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 40 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 41 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 42 43 if (reuse == MAT_REUSE_MATRIX) { 44 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 45 } else { 46 *newmat = B; 47 } 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "MatAXPY_SeqDense" 53 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 54 { 55 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 56 PetscScalar oalpha = alpha; 57 PetscInt j; 58 PetscBLASInt N,m,ldax,lday,one = 1; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 63 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 64 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 65 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 66 if (ldax>m || lday>m) { 67 for (j=0; j<X->cmap->n; j++) { 68 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one)); 69 } 70 } else { 71 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one)); 72 } 73 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 74 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatGetInfo_SeqDense" 80 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 81 { 82 PetscInt N = A->rmap->n*A->cmap->n; 83 84 PetscFunctionBegin; 85 info->block_size = 1.0; 86 info->nz_allocated = (double)N; 87 info->nz_used = (double)N; 88 info->nz_unneeded = (double)0; 89 info->assemblies = (double)A->num_ass; 90 info->mallocs = 0; 91 info->memory = ((PetscObject)A)->mem; 92 info->fill_ratio_given = 0; 93 info->fill_ratio_needed = 0; 94 info->factor_mallocs = 0; 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatScale_SeqDense" 100 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 101 { 102 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 103 PetscScalar oalpha = alpha; 104 PetscErrorCode ierr; 105 PetscBLASInt one = 1,j,nz,lda; 106 107 PetscFunctionBegin; 108 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 109 if (lda>A->rmap->n) { 110 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 111 for (j=0; j<A->cmap->n; j++) { 112 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one)); 113 } 114 } else { 115 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 116 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one)); 117 } 118 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatIsHermitian_SeqDense" 124 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 125 { 126 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 127 PetscInt i,j,m = A->rmap->n,N; 128 PetscScalar *v = a->v; 129 130 PetscFunctionBegin; 131 *fl = PETSC_FALSE; 132 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 133 N = a->lda; 134 135 for (i=0; i<m; i++) { 136 for (j=i+1; j<m; j++) { 137 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 138 } 139 } 140 *fl = PETSC_TRUE; 141 PetscFunctionReturn(0); 142 } 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 146 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 147 { 148 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 149 PetscErrorCode ierr; 150 PetscInt lda = (PetscInt)mat->lda,j,m; 151 152 PetscFunctionBegin; 153 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 154 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 155 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 156 if (cpvalues == MAT_COPY_VALUES) { 157 l = (Mat_SeqDense*)newi->data; 158 if (lda>A->rmap->n) { 159 m = A->rmap->n; 160 for (j=0; j<A->cmap->n; j++) { 161 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 162 } 163 } else { 164 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 165 } 166 } 167 newi->assembled = PETSC_TRUE; 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "MatDuplicate_SeqDense" 173 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 174 { 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 179 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 180 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 181 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 186 extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 187 188 #undef __FUNCT__ 189 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 190 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 191 { 192 MatFactorInfo info; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 197 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 #undef __FUNCT__ 202 #define __FUNCT__ "MatSolve_SeqDense" 203 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 204 { 205 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 206 PetscErrorCode ierr; 207 PetscScalar *x,*y; 208 PetscBLASInt one = 1,info,m; 209 210 PetscFunctionBegin; 211 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 212 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 213 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 214 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 215 if (A->factortype == MAT_FACTOR_LU) { 216 #if defined(PETSC_MISSING_LAPACK_GETRS) 217 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 218 #else 219 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 220 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 221 #endif 222 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 223 #if defined(PETSC_MISSING_LAPACK_POTRS) 224 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 225 #else 226 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 227 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 228 #endif 229 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 230 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 231 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 232 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 #undef __FUNCT__ 237 #define __FUNCT__ "MatMatSolve_SeqDense" 238 PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 239 { 240 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 241 PetscErrorCode ierr; 242 PetscScalar *b,*x; 243 PetscInt n; 244 PetscBLASInt nrhs,info,m; 245 PetscBool flg; 246 247 PetscFunctionBegin; 248 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 249 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 250 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 251 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 252 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 253 254 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 255 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 256 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 257 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 258 259 ierr = PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));CHKERRQ(ierr); 260 261 if (A->factortype == MAT_FACTOR_LU) { 262 #if defined(PETSC_MISSING_LAPACK_GETRS) 263 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 264 #else 265 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 266 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 267 #endif 268 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 269 #if defined(PETSC_MISSING_LAPACK_POTRS) 270 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 271 #else 272 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 273 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 274 #endif 275 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 276 277 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 278 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 279 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "MatSolveTranspose_SeqDense" 285 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 286 { 287 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 288 PetscErrorCode ierr; 289 PetscScalar *x,*y; 290 PetscBLASInt one = 1,info,m; 291 292 PetscFunctionBegin; 293 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 294 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 295 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 296 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 297 /* assume if pivots exist then use LU; else Cholesky */ 298 if (mat->pivots) { 299 #if defined(PETSC_MISSING_LAPACK_GETRS) 300 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 301 #else 302 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 303 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 304 #endif 305 } else { 306 #if defined(PETSC_MISSING_LAPACK_POTRS) 307 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 308 #else 309 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 310 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 311 #endif 312 } 313 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 314 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 315 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatSolveAdd_SeqDense" 321 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 322 { 323 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 324 PetscErrorCode ierr; 325 PetscScalar *x,*y,sone = 1.0; 326 Vec tmp = 0; 327 PetscBLASInt one = 1,info,m; 328 329 PetscFunctionBegin; 330 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 331 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 332 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 333 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 334 if (yy == zz) { 335 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 336 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 337 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 338 } 339 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 340 /* assume if pivots exist then use LU; else Cholesky */ 341 if (mat->pivots) { 342 #if defined(PETSC_MISSING_LAPACK_GETRS) 343 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 344 #else 345 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 346 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 347 #endif 348 } else { 349 #if defined(PETSC_MISSING_LAPACK_POTRS) 350 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 351 #else 352 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 353 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 354 #endif 355 } 356 if (tmp) { 357 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 358 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 359 } else { 360 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 361 } 362 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 363 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 364 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 365 PetscFunctionReturn(0); 366 } 367 368 #undef __FUNCT__ 369 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 370 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 371 { 372 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 373 PetscErrorCode ierr; 374 PetscScalar *x,*y,sone = 1.0; 375 Vec tmp; 376 PetscBLASInt one = 1,info,m; 377 378 PetscFunctionBegin; 379 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 380 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 381 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 382 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 383 if (yy == zz) { 384 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 385 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 386 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 387 } 388 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 389 /* assume if pivots exist then use LU; else Cholesky */ 390 if (mat->pivots) { 391 #if defined(PETSC_MISSING_LAPACK_GETRS) 392 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 393 #else 394 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 395 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 396 #endif 397 } else { 398 #if defined(PETSC_MISSING_LAPACK_POTRS) 399 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 400 #else 401 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 402 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 403 #endif 404 } 405 if (tmp) { 406 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 407 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 408 } else { 409 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 410 } 411 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 412 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 413 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 414 PetscFunctionReturn(0); 415 } 416 417 /* ---------------------------------------------------------------*/ 418 /* COMMENT: I have chosen to hide row permutation in the pivots, 419 rather than put it in the Mat->row slot.*/ 420 #undef __FUNCT__ 421 #define __FUNCT__ "MatLUFactor_SeqDense" 422 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 423 { 424 #if defined(PETSC_MISSING_LAPACK_GETRF) 425 PetscFunctionBegin; 426 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 427 #else 428 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 429 PetscErrorCode ierr; 430 PetscBLASInt n,m,info; 431 432 PetscFunctionBegin; 433 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 434 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 435 if (!mat->pivots) { 436 ierr = PetscMalloc1((A->rmap->n+1),&mat->pivots);CHKERRQ(ierr); 437 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 438 } 439 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 440 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 441 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 442 ierr = PetscFPTrapPop();CHKERRQ(ierr); 443 444 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 445 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 446 A->ops->solve = MatSolve_SeqDense; 447 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 448 A->ops->solveadd = MatSolveAdd_SeqDense; 449 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 450 A->factortype = MAT_FACTOR_LU; 451 452 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 453 #endif 454 PetscFunctionReturn(0); 455 } 456 457 #undef __FUNCT__ 458 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 459 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 460 { 461 #if defined(PETSC_MISSING_LAPACK_POTRF) 462 PetscFunctionBegin; 463 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 464 #else 465 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 466 PetscErrorCode ierr; 467 PetscBLASInt info,n; 468 469 PetscFunctionBegin; 470 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 471 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 472 473 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 474 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 475 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 476 A->ops->solve = MatSolve_SeqDense; 477 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 478 A->ops->solveadd = MatSolveAdd_SeqDense; 479 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 480 A->factortype = MAT_FACTOR_CHOLESKY; 481 482 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 483 #endif 484 PetscFunctionReturn(0); 485 } 486 487 488 #undef __FUNCT__ 489 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 490 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 491 { 492 PetscErrorCode ierr; 493 MatFactorInfo info; 494 495 PetscFunctionBegin; 496 info.fill = 1.0; 497 498 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 499 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 500 PetscFunctionReturn(0); 501 } 502 503 #undef __FUNCT__ 504 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 505 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 506 { 507 PetscFunctionBegin; 508 fact->assembled = PETSC_TRUE; 509 fact->preallocated = PETSC_TRUE; 510 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 511 PetscFunctionReturn(0); 512 } 513 514 #undef __FUNCT__ 515 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 516 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 517 { 518 PetscFunctionBegin; 519 fact->preallocated = PETSC_TRUE; 520 fact->assembled = PETSC_TRUE; 521 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 522 PetscFunctionReturn(0); 523 } 524 525 #undef __FUNCT__ 526 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 527 PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 528 { 529 PetscErrorCode ierr; 530 531 PetscFunctionBegin; 532 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 533 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 534 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 535 if (ftype == MAT_FACTOR_LU) { 536 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 537 } else { 538 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 539 } 540 (*fact)->factortype = ftype; 541 PetscFunctionReturn(0); 542 } 543 544 /* ------------------------------------------------------------------*/ 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatSOR_SeqDense" 547 PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 548 { 549 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 550 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 551 PetscErrorCode ierr; 552 PetscInt m = A->rmap->n,i; 553 PetscBLASInt o = 1,bm; 554 555 PetscFunctionBegin; 556 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 557 if (flag & SOR_ZERO_INITIAL_GUESS) { 558 /* this is a hack fix, should have another version without the second BLASdot */ 559 ierr = VecSet(xx,zero);CHKERRQ(ierr); 560 } 561 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 562 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 563 its = its*lits; 564 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 565 while (its--) { 566 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 567 for (i=0; i<m; i++) { 568 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 569 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 570 } 571 } 572 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 573 for (i=m-1; i>=0; i--) { 574 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 575 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 576 } 577 } 578 } 579 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 580 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 /* -----------------------------------------------------------------*/ 585 #undef __FUNCT__ 586 #define __FUNCT__ "MatMultTranspose_SeqDense" 587 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 588 { 589 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 590 PetscScalar *v = mat->v,*x,*y; 591 PetscErrorCode ierr; 592 PetscBLASInt m, n,_One=1; 593 PetscScalar _DOne=1.0,_DZero=0.0; 594 595 PetscFunctionBegin; 596 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 597 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 598 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 599 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 600 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 601 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 602 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 603 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 604 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 605 PetscFunctionReturn(0); 606 } 607 608 #undef __FUNCT__ 609 #define __FUNCT__ "MatMult_SeqDense" 610 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 611 { 612 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 613 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 614 PetscErrorCode ierr; 615 PetscBLASInt m, n, _One=1; 616 617 PetscFunctionBegin; 618 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 619 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 620 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 621 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 622 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 623 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 624 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 625 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 626 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 627 PetscFunctionReturn(0); 628 } 629 630 #undef __FUNCT__ 631 #define __FUNCT__ "MatMultAdd_SeqDense" 632 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 633 { 634 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 635 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 636 PetscErrorCode ierr; 637 PetscBLASInt m, n, _One=1; 638 639 PetscFunctionBegin; 640 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 641 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 642 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 643 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 644 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 645 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 646 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 647 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 648 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 649 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 650 PetscFunctionReturn(0); 651 } 652 653 #undef __FUNCT__ 654 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 655 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 656 { 657 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 658 PetscScalar *v = mat->v,*x,*y; 659 PetscErrorCode ierr; 660 PetscBLASInt m, n, _One=1; 661 PetscScalar _DOne=1.0; 662 663 PetscFunctionBegin; 664 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 665 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 666 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 667 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 668 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 669 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 670 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 671 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 672 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 673 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 674 PetscFunctionReturn(0); 675 } 676 677 /* -----------------------------------------------------------------*/ 678 #undef __FUNCT__ 679 #define __FUNCT__ "MatGetRow_SeqDense" 680 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 681 { 682 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 683 PetscScalar *v; 684 PetscErrorCode ierr; 685 PetscInt i; 686 687 PetscFunctionBegin; 688 *ncols = A->cmap->n; 689 if (cols) { 690 ierr = PetscMalloc1((A->cmap->n+1),cols);CHKERRQ(ierr); 691 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 692 } 693 if (vals) { 694 ierr = PetscMalloc1((A->cmap->n+1),vals);CHKERRQ(ierr); 695 v = mat->v + row; 696 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 697 } 698 PetscFunctionReturn(0); 699 } 700 701 #undef __FUNCT__ 702 #define __FUNCT__ "MatRestoreRow_SeqDense" 703 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 704 { 705 PetscErrorCode ierr; 706 707 PetscFunctionBegin; 708 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 709 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 710 PetscFunctionReturn(0); 711 } 712 /* ----------------------------------------------------------------*/ 713 #undef __FUNCT__ 714 #define __FUNCT__ "MatSetValues_SeqDense" 715 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 716 { 717 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 718 PetscInt i,j,idx=0; 719 720 PetscFunctionBegin; 721 if (!mat->roworiented) { 722 if (addv == INSERT_VALUES) { 723 for (j=0; j<n; j++) { 724 if (indexn[j] < 0) {idx += m; continue;} 725 #if defined(PETSC_USE_DEBUG) 726 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 727 #endif 728 for (i=0; i<m; i++) { 729 if (indexm[i] < 0) {idx++; continue;} 730 #if defined(PETSC_USE_DEBUG) 731 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 732 #endif 733 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 734 } 735 } 736 } else { 737 for (j=0; j<n; j++) { 738 if (indexn[j] < 0) {idx += m; continue;} 739 #if defined(PETSC_USE_DEBUG) 740 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 741 #endif 742 for (i=0; i<m; i++) { 743 if (indexm[i] < 0) {idx++; continue;} 744 #if defined(PETSC_USE_DEBUG) 745 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 746 #endif 747 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 748 } 749 } 750 } 751 } else { 752 if (addv == INSERT_VALUES) { 753 for (i=0; i<m; i++) { 754 if (indexm[i] < 0) { idx += n; continue;} 755 #if defined(PETSC_USE_DEBUG) 756 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 757 #endif 758 for (j=0; j<n; j++) { 759 if (indexn[j] < 0) { idx++; continue;} 760 #if defined(PETSC_USE_DEBUG) 761 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 762 #endif 763 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 764 } 765 } 766 } else { 767 for (i=0; i<m; i++) { 768 if (indexm[i] < 0) { idx += n; continue;} 769 #if defined(PETSC_USE_DEBUG) 770 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 771 #endif 772 for (j=0; j<n; j++) { 773 if (indexn[j] < 0) { idx++; continue;} 774 #if defined(PETSC_USE_DEBUG) 775 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 776 #endif 777 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 778 } 779 } 780 } 781 } 782 PetscFunctionReturn(0); 783 } 784 785 #undef __FUNCT__ 786 #define __FUNCT__ "MatGetValues_SeqDense" 787 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 788 { 789 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 790 PetscInt i,j; 791 792 PetscFunctionBegin; 793 /* row-oriented output */ 794 for (i=0; i<m; i++) { 795 if (indexm[i] < 0) {v += n;continue;} 796 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 797 for (j=0; j<n; j++) { 798 if (indexn[j] < 0) {v++; continue;} 799 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 800 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 801 } 802 } 803 PetscFunctionReturn(0); 804 } 805 806 /* -----------------------------------------------------------------*/ 807 808 #undef __FUNCT__ 809 #define __FUNCT__ "MatLoad_SeqDense" 810 PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 811 { 812 Mat_SeqDense *a; 813 PetscErrorCode ierr; 814 PetscInt *scols,i,j,nz,header[4]; 815 int fd; 816 PetscMPIInt size; 817 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 818 PetscScalar *vals,*svals,*v,*w; 819 MPI_Comm comm; 820 821 PetscFunctionBegin; 822 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 823 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 824 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 825 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 826 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 827 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 828 M = header[1]; N = header[2]; nz = header[3]; 829 830 /* set global size if not set already*/ 831 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 832 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 833 } else { 834 /* if sizes and type are already set, check if the vector global sizes are correct */ 835 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 836 if (M != grows || N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols); 837 } 838 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 839 840 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 841 a = (Mat_SeqDense*)newmat->data; 842 v = a->v; 843 /* Allocate some temp space to read in the values and then flip them 844 from row major to column major */ 845 ierr = PetscMalloc1((M*N > 0 ? M*N : 1),&w);CHKERRQ(ierr); 846 /* read in nonzero values */ 847 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 848 /* now flip the values and store them in the matrix*/ 849 for (j=0; j<N; j++) { 850 for (i=0; i<M; i++) { 851 *v++ =w[i*N+j]; 852 } 853 } 854 ierr = PetscFree(w);CHKERRQ(ierr); 855 } else { 856 /* read row lengths */ 857 ierr = PetscMalloc1((M+1),&rowlengths);CHKERRQ(ierr); 858 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 859 860 a = (Mat_SeqDense*)newmat->data; 861 v = a->v; 862 863 /* read column indices and nonzeros */ 864 ierr = PetscMalloc1((nz+1),&scols);CHKERRQ(ierr); 865 cols = scols; 866 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 867 ierr = PetscMalloc1((nz+1),&svals);CHKERRQ(ierr); 868 vals = svals; 869 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 870 871 /* insert into matrix */ 872 for (i=0; i<M; i++) { 873 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 874 svals += rowlengths[i]; scols += rowlengths[i]; 875 } 876 ierr = PetscFree(vals);CHKERRQ(ierr); 877 ierr = PetscFree(cols);CHKERRQ(ierr); 878 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 879 } 880 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 881 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 882 PetscFunctionReturn(0); 883 } 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "MatView_SeqDense_ASCII" 887 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 888 { 889 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 890 PetscErrorCode ierr; 891 PetscInt i,j; 892 const char *name; 893 PetscScalar *v; 894 PetscViewerFormat format; 895 #if defined(PETSC_USE_COMPLEX) 896 PetscBool allreal = PETSC_TRUE; 897 #endif 898 899 PetscFunctionBegin; 900 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 901 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 902 PetscFunctionReturn(0); /* do nothing for now */ 903 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 904 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 905 for (i=0; i<A->rmap->n; i++) { 906 v = a->v + i; 907 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 908 for (j=0; j<A->cmap->n; j++) { 909 #if defined(PETSC_USE_COMPLEX) 910 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 911 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 912 } else if (PetscRealPart(*v)) { 913 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 914 } 915 #else 916 if (*v) { 917 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 918 } 919 #endif 920 v += a->lda; 921 } 922 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 923 } 924 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 925 } else { 926 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 927 #if defined(PETSC_USE_COMPLEX) 928 /* determine if matrix has all real values */ 929 v = a->v; 930 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 931 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 932 } 933 #endif 934 if (format == PETSC_VIEWER_ASCII_MATLAB) { 935 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 936 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 937 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 938 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 939 } 940 941 for (i=0; i<A->rmap->n; i++) { 942 v = a->v + i; 943 for (j=0; j<A->cmap->n; j++) { 944 #if defined(PETSC_USE_COMPLEX) 945 if (allreal) { 946 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 947 } else { 948 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 949 } 950 #else 951 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 952 #endif 953 v += a->lda; 954 } 955 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 956 } 957 if (format == PETSC_VIEWER_ASCII_MATLAB) { 958 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 959 } 960 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 961 } 962 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 963 PetscFunctionReturn(0); 964 } 965 966 #undef __FUNCT__ 967 #define __FUNCT__ "MatView_SeqDense_Binary" 968 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 969 { 970 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 971 PetscErrorCode ierr; 972 int fd; 973 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 974 PetscScalar *v,*anonz,*vals; 975 PetscViewerFormat format; 976 977 PetscFunctionBegin; 978 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 979 980 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 981 if (format == PETSC_VIEWER_NATIVE) { 982 /* store the matrix as a dense matrix */ 983 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 984 985 col_lens[0] = MAT_FILE_CLASSID; 986 col_lens[1] = m; 987 col_lens[2] = n; 988 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 989 990 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 991 ierr = PetscFree(col_lens);CHKERRQ(ierr); 992 993 /* write out matrix, by rows */ 994 ierr = PetscMalloc1((m*n+1),&vals);CHKERRQ(ierr); 995 v = a->v; 996 for (j=0; j<n; j++) { 997 for (i=0; i<m; i++) { 998 vals[j + i*n] = *v++; 999 } 1000 } 1001 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1002 ierr = PetscFree(vals);CHKERRQ(ierr); 1003 } else { 1004 ierr = PetscMalloc1((4+nz),&col_lens);CHKERRQ(ierr); 1005 1006 col_lens[0] = MAT_FILE_CLASSID; 1007 col_lens[1] = m; 1008 col_lens[2] = n; 1009 col_lens[3] = nz; 1010 1011 /* store lengths of each row and write (including header) to file */ 1012 for (i=0; i<m; i++) col_lens[4+i] = n; 1013 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1014 1015 /* Possibly should write in smaller increments, not whole matrix at once? */ 1016 /* store column indices (zero start index) */ 1017 ict = 0; 1018 for (i=0; i<m; i++) { 1019 for (j=0; j<n; j++) col_lens[ict++] = j; 1020 } 1021 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1022 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1023 1024 /* store nonzero values */ 1025 ierr = PetscMalloc1((nz+1),&anonz);CHKERRQ(ierr); 1026 ict = 0; 1027 for (i=0; i<m; i++) { 1028 v = a->v + i; 1029 for (j=0; j<n; j++) { 1030 anonz[ict++] = *v; v += a->lda; 1031 } 1032 } 1033 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1034 ierr = PetscFree(anonz);CHKERRQ(ierr); 1035 } 1036 PetscFunctionReturn(0); 1037 } 1038 1039 #include <petscdraw.h> 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1042 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1043 { 1044 Mat A = (Mat) Aa; 1045 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1046 PetscErrorCode ierr; 1047 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1048 PetscScalar *v = a->v; 1049 PetscViewer viewer; 1050 PetscDraw popup; 1051 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1052 PetscViewerFormat format; 1053 1054 PetscFunctionBegin; 1055 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1056 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1057 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1058 1059 /* Loop over matrix elements drawing boxes */ 1060 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1061 /* Blue for negative and Red for positive */ 1062 color = PETSC_DRAW_BLUE; 1063 for (j = 0; j < n; j++) { 1064 x_l = j; 1065 x_r = x_l + 1.0; 1066 for (i = 0; i < m; i++) { 1067 y_l = m - i - 1.0; 1068 y_r = y_l + 1.0; 1069 if (PetscRealPart(v[j*m+i]) > 0.) { 1070 color = PETSC_DRAW_RED; 1071 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1072 color = PETSC_DRAW_BLUE; 1073 } else { 1074 continue; 1075 } 1076 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1077 } 1078 } 1079 } else { 1080 /* use contour shading to indicate magnitude of values */ 1081 /* first determine max of all nonzero values */ 1082 for (i = 0; i < m*n; i++) { 1083 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1084 } 1085 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1086 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1087 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1088 for (j = 0; j < n; j++) { 1089 x_l = j; 1090 x_r = x_l + 1.0; 1091 for (i = 0; i < m; i++) { 1092 y_l = m - i - 1.0; 1093 y_r = y_l + 1.0; 1094 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1095 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1096 } 1097 } 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNCT__ 1103 #define __FUNCT__ "MatView_SeqDense_Draw" 1104 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1105 { 1106 PetscDraw draw; 1107 PetscBool isnull; 1108 PetscReal xr,yr,xl,yl,h,w; 1109 PetscErrorCode ierr; 1110 1111 PetscFunctionBegin; 1112 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1113 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1114 if (isnull) PetscFunctionReturn(0); 1115 1116 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1117 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1118 xr += w; yr += h; xl = -w; yl = -h; 1119 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1120 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1121 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "MatView_SeqDense" 1127 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1128 { 1129 PetscErrorCode ierr; 1130 PetscBool iascii,isbinary,isdraw; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1134 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1135 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1136 1137 if (iascii) { 1138 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1139 } else if (isbinary) { 1140 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1141 } else if (isdraw) { 1142 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1143 } 1144 PetscFunctionReturn(0); 1145 } 1146 1147 #undef __FUNCT__ 1148 #define __FUNCT__ "MatDestroy_SeqDense" 1149 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1150 { 1151 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1152 PetscErrorCode ierr; 1153 1154 PetscFunctionBegin; 1155 #if defined(PETSC_USE_LOG) 1156 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1157 #endif 1158 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1159 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1160 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1161 1162 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1163 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1164 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1171 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1172 PetscFunctionReturn(0); 1173 } 1174 1175 #undef __FUNCT__ 1176 #define __FUNCT__ "MatTranspose_SeqDense" 1177 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1178 { 1179 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1180 PetscErrorCode ierr; 1181 PetscInt k,j,m,n,M; 1182 PetscScalar *v,tmp; 1183 1184 PetscFunctionBegin; 1185 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1186 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1187 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1188 else { 1189 for (j=0; j<m; j++) { 1190 for (k=0; k<j; k++) { 1191 tmp = v[j + k*M]; 1192 v[j + k*M] = v[k + j*M]; 1193 v[k + j*M] = tmp; 1194 } 1195 } 1196 } 1197 } else { /* out-of-place transpose */ 1198 Mat tmat; 1199 Mat_SeqDense *tmatd; 1200 PetscScalar *v2; 1201 1202 if (reuse == MAT_INITIAL_MATRIX) { 1203 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1204 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1205 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1206 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1207 } else { 1208 tmat = *matout; 1209 } 1210 tmatd = (Mat_SeqDense*)tmat->data; 1211 v = mat->v; v2 = tmatd->v; 1212 for (j=0; j<n; j++) { 1213 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1214 } 1215 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1216 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1217 *matout = tmat; 1218 } 1219 PetscFunctionReturn(0); 1220 } 1221 1222 #undef __FUNCT__ 1223 #define __FUNCT__ "MatEqual_SeqDense" 1224 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1225 { 1226 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1227 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1228 PetscInt i,j; 1229 PetscScalar *v1,*v2; 1230 1231 PetscFunctionBegin; 1232 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1233 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1234 for (i=0; i<A1->rmap->n; i++) { 1235 v1 = mat1->v+i; v2 = mat2->v+i; 1236 for (j=0; j<A1->cmap->n; j++) { 1237 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1238 v1 += mat1->lda; v2 += mat2->lda; 1239 } 1240 } 1241 *flg = PETSC_TRUE; 1242 PetscFunctionReturn(0); 1243 } 1244 1245 #undef __FUNCT__ 1246 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1247 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1248 { 1249 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1250 PetscErrorCode ierr; 1251 PetscInt i,n,len; 1252 PetscScalar *x,zero = 0.0; 1253 1254 PetscFunctionBegin; 1255 ierr = VecSet(v,zero);CHKERRQ(ierr); 1256 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1257 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1258 len = PetscMin(A->rmap->n,A->cmap->n); 1259 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1260 for (i=0; i<len; i++) { 1261 x[i] = mat->v[i*mat->lda + i]; 1262 } 1263 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 #undef __FUNCT__ 1268 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1269 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1270 { 1271 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1272 PetscScalar *l,*r,x,*v; 1273 PetscErrorCode ierr; 1274 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1275 1276 PetscFunctionBegin; 1277 if (ll) { 1278 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1279 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1280 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1281 for (i=0; i<m; i++) { 1282 x = l[i]; 1283 v = mat->v + i; 1284 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1285 } 1286 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1287 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1288 } 1289 if (rr) { 1290 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1291 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1292 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1293 for (i=0; i<n; i++) { 1294 x = r[i]; 1295 v = mat->v + i*m; 1296 for (j=0; j<m; j++) (*v++) *= x; 1297 } 1298 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1299 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1300 } 1301 PetscFunctionReturn(0); 1302 } 1303 1304 #undef __FUNCT__ 1305 #define __FUNCT__ "MatNorm_SeqDense" 1306 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1307 { 1308 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1309 PetscScalar *v = mat->v; 1310 PetscReal sum = 0.0; 1311 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1312 PetscErrorCode ierr; 1313 1314 PetscFunctionBegin; 1315 if (type == NORM_FROBENIUS) { 1316 if (lda>m) { 1317 for (j=0; j<A->cmap->n; j++) { 1318 v = mat->v+j*lda; 1319 for (i=0; i<m; i++) { 1320 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1321 } 1322 } 1323 } else { 1324 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1325 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1326 } 1327 } 1328 *nrm = PetscSqrtReal(sum); 1329 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1330 } else if (type == NORM_1) { 1331 *nrm = 0.0; 1332 for (j=0; j<A->cmap->n; j++) { 1333 v = mat->v + j*mat->lda; 1334 sum = 0.0; 1335 for (i=0; i<A->rmap->n; i++) { 1336 sum += PetscAbsScalar(*v); v++; 1337 } 1338 if (sum > *nrm) *nrm = sum; 1339 } 1340 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1341 } else if (type == NORM_INFINITY) { 1342 *nrm = 0.0; 1343 for (j=0; j<A->rmap->n; j++) { 1344 v = mat->v + j; 1345 sum = 0.0; 1346 for (i=0; i<A->cmap->n; i++) { 1347 sum += PetscAbsScalar(*v); v += mat->lda; 1348 } 1349 if (sum > *nrm) *nrm = sum; 1350 } 1351 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1352 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1353 PetscFunctionReturn(0); 1354 } 1355 1356 #undef __FUNCT__ 1357 #define __FUNCT__ "MatSetOption_SeqDense" 1358 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1359 { 1360 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1361 PetscErrorCode ierr; 1362 1363 PetscFunctionBegin; 1364 switch (op) { 1365 case MAT_ROW_ORIENTED: 1366 aij->roworiented = flg; 1367 break; 1368 case MAT_NEW_NONZERO_LOCATIONS: 1369 case MAT_NEW_NONZERO_LOCATION_ERR: 1370 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1371 case MAT_NEW_DIAGONALS: 1372 case MAT_KEEP_NONZERO_PATTERN: 1373 case MAT_IGNORE_OFF_PROC_ENTRIES: 1374 case MAT_USE_HASH_TABLE: 1375 case MAT_IGNORE_LOWER_TRIANGULAR: 1376 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1377 break; 1378 case MAT_SPD: 1379 case MAT_SYMMETRIC: 1380 case MAT_STRUCTURALLY_SYMMETRIC: 1381 case MAT_HERMITIAN: 1382 case MAT_SYMMETRY_ETERNAL: 1383 /* These options are handled directly by MatSetOption() */ 1384 break; 1385 default: 1386 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1387 } 1388 PetscFunctionReturn(0); 1389 } 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "MatZeroEntries_SeqDense" 1393 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1394 { 1395 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1396 PetscErrorCode ierr; 1397 PetscInt lda=l->lda,m=A->rmap->n,j; 1398 1399 PetscFunctionBegin; 1400 if (lda>m) { 1401 for (j=0; j<A->cmap->n; j++) { 1402 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1403 } 1404 } else { 1405 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1406 } 1407 PetscFunctionReturn(0); 1408 } 1409 1410 #undef __FUNCT__ 1411 #define __FUNCT__ "MatZeroRows_SeqDense" 1412 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1413 { 1414 PetscErrorCode ierr; 1415 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1416 PetscInt m = l->lda, n = A->cmap->n, i,j; 1417 PetscScalar *slot,*bb; 1418 const PetscScalar *xx; 1419 1420 PetscFunctionBegin; 1421 #if defined(PETSC_USE_DEBUG) 1422 for (i=0; i<N; i++) { 1423 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1424 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 1425 } 1426 #endif 1427 1428 /* fix right hand side if needed */ 1429 if (x && b) { 1430 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1431 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1432 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1433 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1434 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1435 } 1436 1437 for (i=0; i<N; i++) { 1438 slot = l->v + rows[i]; 1439 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1440 } 1441 if (diag != 0.0) { 1442 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1443 for (i=0; i<N; i++) { 1444 slot = l->v + (m+1)*rows[i]; 1445 *slot = diag; 1446 } 1447 } 1448 PetscFunctionReturn(0); 1449 } 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1453 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1454 { 1455 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1456 1457 PetscFunctionBegin; 1458 if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1459 *array = mat->v; 1460 PetscFunctionReturn(0); 1461 } 1462 1463 #undef __FUNCT__ 1464 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1465 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1466 { 1467 PetscFunctionBegin; 1468 *array = 0; /* user cannot accidently use the array later */ 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "MatDenseGetArray" 1474 /*@C 1475 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1476 1477 Not Collective 1478 1479 Input Parameter: 1480 . mat - a MATSEQDENSE matrix 1481 1482 Output Parameter: 1483 . array - pointer to the data 1484 1485 Level: intermediate 1486 1487 .seealso: MatDenseRestoreArray() 1488 @*/ 1489 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1490 { 1491 PetscErrorCode ierr; 1492 1493 PetscFunctionBegin; 1494 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 #undef __FUNCT__ 1499 #define __FUNCT__ "MatDenseRestoreArray" 1500 /*@C 1501 MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray() 1502 1503 Not Collective 1504 1505 Input Parameters: 1506 . mat - a MATSEQDENSE matrix 1507 . array - pointer to the data 1508 1509 Level: intermediate 1510 1511 .seealso: MatDenseGetArray() 1512 @*/ 1513 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1514 { 1515 PetscErrorCode ierr; 1516 1517 PetscFunctionBegin; 1518 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1519 PetscFunctionReturn(0); 1520 } 1521 1522 #undef __FUNCT__ 1523 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1524 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1525 { 1526 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1527 PetscErrorCode ierr; 1528 PetscInt i,j,nrows,ncols; 1529 const PetscInt *irow,*icol; 1530 PetscScalar *av,*bv,*v = mat->v; 1531 Mat newmat; 1532 1533 PetscFunctionBegin; 1534 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1535 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1536 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1537 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1538 1539 /* Check submatrixcall */ 1540 if (scall == MAT_REUSE_MATRIX) { 1541 PetscInt n_cols,n_rows; 1542 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1543 if (n_rows != nrows || n_cols != ncols) { 1544 /* resize the result matrix to match number of requested rows/columns */ 1545 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1546 } 1547 newmat = *B; 1548 } else { 1549 /* Create and fill new matrix */ 1550 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1551 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1552 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1553 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1554 } 1555 1556 /* Now extract the data pointers and do the copy,column at a time */ 1557 bv = ((Mat_SeqDense*)newmat->data)->v; 1558 1559 for (i=0; i<ncols; i++) { 1560 av = v + mat->lda*icol[i]; 1561 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1562 } 1563 1564 /* Assemble the matrices so that the correct flags are set */ 1565 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1566 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1567 1568 /* Free work space */ 1569 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1570 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1571 *B = newmat; 1572 PetscFunctionReturn(0); 1573 } 1574 1575 #undef __FUNCT__ 1576 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1577 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1578 { 1579 PetscErrorCode ierr; 1580 PetscInt i; 1581 1582 PetscFunctionBegin; 1583 if (scall == MAT_INITIAL_MATRIX) { 1584 ierr = PetscMalloc1((n+1),B);CHKERRQ(ierr); 1585 } 1586 1587 for (i=0; i<n; i++) { 1588 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1589 } 1590 PetscFunctionReturn(0); 1591 } 1592 1593 #undef __FUNCT__ 1594 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1595 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1596 { 1597 PetscFunctionBegin; 1598 PetscFunctionReturn(0); 1599 } 1600 1601 #undef __FUNCT__ 1602 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1603 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1604 { 1605 PetscFunctionBegin; 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "MatCopy_SeqDense" 1611 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1612 { 1613 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1614 PetscErrorCode ierr; 1615 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1616 1617 PetscFunctionBegin; 1618 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1619 if (A->ops->copy != B->ops->copy) { 1620 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1621 PetscFunctionReturn(0); 1622 } 1623 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1624 if (lda1>m || lda2>m) { 1625 for (j=0; j<n; j++) { 1626 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1627 } 1628 } else { 1629 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1630 } 1631 PetscFunctionReturn(0); 1632 } 1633 1634 #undef __FUNCT__ 1635 #define __FUNCT__ "MatSetUp_SeqDense" 1636 PetscErrorCode MatSetUp_SeqDense(Mat A) 1637 { 1638 PetscErrorCode ierr; 1639 1640 PetscFunctionBegin; 1641 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 #undef __FUNCT__ 1646 #define __FUNCT__ "MatConjugate_SeqDense" 1647 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1648 { 1649 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1650 PetscInt i,nz = A->rmap->n*A->cmap->n; 1651 PetscScalar *aa = a->v; 1652 1653 PetscFunctionBegin; 1654 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1655 PetscFunctionReturn(0); 1656 } 1657 1658 #undef __FUNCT__ 1659 #define __FUNCT__ "MatRealPart_SeqDense" 1660 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1661 { 1662 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1663 PetscInt i,nz = A->rmap->n*A->cmap->n; 1664 PetscScalar *aa = a->v; 1665 1666 PetscFunctionBegin; 1667 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1673 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1674 { 1675 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1676 PetscInt i,nz = A->rmap->n*A->cmap->n; 1677 PetscScalar *aa = a->v; 1678 1679 PetscFunctionBegin; 1680 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 /* ----------------------------------------------------------------*/ 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1687 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1688 { 1689 PetscErrorCode ierr; 1690 1691 PetscFunctionBegin; 1692 if (scall == MAT_INITIAL_MATRIX) { 1693 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1694 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1695 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1696 } 1697 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1698 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1699 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1705 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1706 { 1707 PetscErrorCode ierr; 1708 PetscInt m=A->rmap->n,n=B->cmap->n; 1709 Mat Cmat; 1710 1711 PetscFunctionBegin; 1712 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 1713 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1714 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1715 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1716 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1717 1718 *C = Cmat; 1719 PetscFunctionReturn(0); 1720 } 1721 1722 #undef __FUNCT__ 1723 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1724 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1725 { 1726 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1727 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1728 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1729 PetscBLASInt m,n,k; 1730 PetscScalar _DOne=1.0,_DZero=0.0; 1731 PetscErrorCode ierr; 1732 1733 PetscFunctionBegin; 1734 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1735 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1736 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1737 PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1743 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1744 { 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 if (scall == MAT_INITIAL_MATRIX) { 1749 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1750 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1751 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1752 } 1753 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1754 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1755 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 #undef __FUNCT__ 1760 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1761 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1762 { 1763 PetscErrorCode ierr; 1764 PetscInt m=A->cmap->n,n=B->cmap->n; 1765 Mat Cmat; 1766 1767 PetscFunctionBegin; 1768 if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n); 1769 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1770 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1771 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1772 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1773 1774 Cmat->assembled = PETSC_TRUE; 1775 1776 *C = Cmat; 1777 PetscFunctionReturn(0); 1778 } 1779 1780 #undef __FUNCT__ 1781 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1782 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1783 { 1784 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1785 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1786 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1787 PetscBLASInt m,n,k; 1788 PetscScalar _DOne=1.0,_DZero=0.0; 1789 PetscErrorCode ierr; 1790 1791 PetscFunctionBegin; 1792 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1793 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1794 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1795 /* 1796 Note the m and n arguments below are the number rows and columns of A', not A! 1797 */ 1798 PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1799 PetscFunctionReturn(0); 1800 } 1801 1802 #undef __FUNCT__ 1803 #define __FUNCT__ "MatGetRowMax_SeqDense" 1804 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1805 { 1806 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1807 PetscErrorCode ierr; 1808 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1809 PetscScalar *x; 1810 MatScalar *aa = a->v; 1811 1812 PetscFunctionBegin; 1813 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1814 1815 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1816 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1817 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1818 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1819 for (i=0; i<m; i++) { 1820 x[i] = aa[i]; if (idx) idx[i] = 0; 1821 for (j=1; j<n; j++) { 1822 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1823 } 1824 } 1825 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1831 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1832 { 1833 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1834 PetscErrorCode ierr; 1835 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1836 PetscScalar *x; 1837 PetscReal atmp; 1838 MatScalar *aa = a->v; 1839 1840 PetscFunctionBegin; 1841 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1842 1843 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1844 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1845 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1846 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1847 for (i=0; i<m; i++) { 1848 x[i] = PetscAbsScalar(aa[i]); 1849 for (j=1; j<n; j++) { 1850 atmp = PetscAbsScalar(aa[i+m*j]); 1851 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1852 } 1853 } 1854 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1855 PetscFunctionReturn(0); 1856 } 1857 1858 #undef __FUNCT__ 1859 #define __FUNCT__ "MatGetRowMin_SeqDense" 1860 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1861 { 1862 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1863 PetscErrorCode ierr; 1864 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1865 PetscScalar *x; 1866 MatScalar *aa = a->v; 1867 1868 PetscFunctionBegin; 1869 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1870 1871 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1872 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1873 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1874 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1875 for (i=0; i<m; i++) { 1876 x[i] = aa[i]; if (idx) idx[i] = 0; 1877 for (j=1; j<n; j++) { 1878 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1879 } 1880 } 1881 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1882 PetscFunctionReturn(0); 1883 } 1884 1885 #undef __FUNCT__ 1886 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1887 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1888 { 1889 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1890 PetscErrorCode ierr; 1891 PetscScalar *x; 1892 1893 PetscFunctionBegin; 1894 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1895 1896 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1897 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1898 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 1903 #undef __FUNCT__ 1904 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1905 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1906 { 1907 PetscErrorCode ierr; 1908 PetscInt i,j,m,n; 1909 PetscScalar *a; 1910 1911 PetscFunctionBegin; 1912 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1913 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1914 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1915 if (type == NORM_2) { 1916 for (i=0; i<n; i++) { 1917 for (j=0; j<m; j++) { 1918 norms[i] += PetscAbsScalar(a[j]*a[j]); 1919 } 1920 a += m; 1921 } 1922 } else if (type == NORM_1) { 1923 for (i=0; i<n; i++) { 1924 for (j=0; j<m; j++) { 1925 norms[i] += PetscAbsScalar(a[j]); 1926 } 1927 a += m; 1928 } 1929 } else if (type == NORM_INFINITY) { 1930 for (i=0; i<n; i++) { 1931 for (j=0; j<m; j++) { 1932 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1933 } 1934 a += m; 1935 } 1936 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1937 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1938 if (type == NORM_2) { 1939 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1940 } 1941 PetscFunctionReturn(0); 1942 } 1943 1944 #undef __FUNCT__ 1945 #define __FUNCT__ "MatSetRandom_SeqDense" 1946 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1947 { 1948 PetscErrorCode ierr; 1949 PetscScalar *a; 1950 PetscInt m,n,i; 1951 1952 PetscFunctionBegin; 1953 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1954 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1955 for (i=0; i<m*n; i++) { 1956 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1957 } 1958 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1959 PetscFunctionReturn(0); 1960 } 1961 1962 1963 /* -------------------------------------------------------------------*/ 1964 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 1965 MatGetRow_SeqDense, 1966 MatRestoreRow_SeqDense, 1967 MatMult_SeqDense, 1968 /* 4*/ MatMultAdd_SeqDense, 1969 MatMultTranspose_SeqDense, 1970 MatMultTransposeAdd_SeqDense, 1971 0, 1972 0, 1973 0, 1974 /* 10*/ 0, 1975 MatLUFactor_SeqDense, 1976 MatCholeskyFactor_SeqDense, 1977 MatSOR_SeqDense, 1978 MatTranspose_SeqDense, 1979 /* 15*/ MatGetInfo_SeqDense, 1980 MatEqual_SeqDense, 1981 MatGetDiagonal_SeqDense, 1982 MatDiagonalScale_SeqDense, 1983 MatNorm_SeqDense, 1984 /* 20*/ MatAssemblyBegin_SeqDense, 1985 MatAssemblyEnd_SeqDense, 1986 MatSetOption_SeqDense, 1987 MatZeroEntries_SeqDense, 1988 /* 24*/ MatZeroRows_SeqDense, 1989 0, 1990 0, 1991 0, 1992 0, 1993 /* 29*/ MatSetUp_SeqDense, 1994 0, 1995 0, 1996 0, 1997 0, 1998 /* 34*/ MatDuplicate_SeqDense, 1999 0, 2000 0, 2001 0, 2002 0, 2003 /* 39*/ MatAXPY_SeqDense, 2004 MatGetSubMatrices_SeqDense, 2005 0, 2006 MatGetValues_SeqDense, 2007 MatCopy_SeqDense, 2008 /* 44*/ MatGetRowMax_SeqDense, 2009 MatScale_SeqDense, 2010 0, 2011 0, 2012 0, 2013 /* 49*/ MatSetRandom_SeqDense, 2014 0, 2015 0, 2016 0, 2017 0, 2018 /* 54*/ 0, 2019 0, 2020 0, 2021 0, 2022 0, 2023 /* 59*/ 0, 2024 MatDestroy_SeqDense, 2025 MatView_SeqDense, 2026 0, 2027 0, 2028 /* 64*/ 0, 2029 0, 2030 0, 2031 0, 2032 0, 2033 /* 69*/ MatGetRowMaxAbs_SeqDense, 2034 0, 2035 0, 2036 0, 2037 0, 2038 /* 74*/ 0, 2039 0, 2040 0, 2041 0, 2042 0, 2043 /* 79*/ 0, 2044 0, 2045 0, 2046 0, 2047 /* 83*/ MatLoad_SeqDense, 2048 0, 2049 MatIsHermitian_SeqDense, 2050 0, 2051 0, 2052 0, 2053 /* 89*/ MatMatMult_SeqDense_SeqDense, 2054 MatMatMultSymbolic_SeqDense_SeqDense, 2055 MatMatMultNumeric_SeqDense_SeqDense, 2056 0, 2057 0, 2058 /* 94*/ 0, 2059 0, 2060 0, 2061 0, 2062 0, 2063 /* 99*/ 0, 2064 0, 2065 0, 2066 MatConjugate_SeqDense, 2067 0, 2068 /*104*/ 0, 2069 MatRealPart_SeqDense, 2070 MatImaginaryPart_SeqDense, 2071 0, 2072 0, 2073 /*109*/ MatMatSolve_SeqDense, 2074 0, 2075 MatGetRowMin_SeqDense, 2076 MatGetColumnVector_SeqDense, 2077 0, 2078 /*114*/ 0, 2079 0, 2080 0, 2081 0, 2082 0, 2083 /*119*/ 0, 2084 0, 2085 0, 2086 0, 2087 0, 2088 /*124*/ 0, 2089 MatGetColumnNorms_SeqDense, 2090 0, 2091 0, 2092 0, 2093 /*129*/ 0, 2094 MatTransposeMatMult_SeqDense_SeqDense, 2095 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2096 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2097 0, 2098 /*134*/ 0, 2099 0, 2100 0, 2101 0, 2102 0, 2103 /*139*/ 0, 2104 0, 2105 0 2106 }; 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "MatCreateSeqDense" 2110 /*@C 2111 MatCreateSeqDense - Creates a sequential dense matrix that 2112 is stored in column major order (the usual Fortran 77 manner). Many 2113 of the matrix operations use the BLAS and LAPACK routines. 2114 2115 Collective on MPI_Comm 2116 2117 Input Parameters: 2118 + comm - MPI communicator, set to PETSC_COMM_SELF 2119 . m - number of rows 2120 . n - number of columns 2121 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2122 to control all matrix memory allocation. 2123 2124 Output Parameter: 2125 . A - the matrix 2126 2127 Notes: 2128 The data input variable is intended primarily for Fortran programmers 2129 who wish to allocate their own matrix memory space. Most users should 2130 set data=NULL. 2131 2132 Level: intermediate 2133 2134 .keywords: dense, matrix, LAPACK, BLAS 2135 2136 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2137 @*/ 2138 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2139 { 2140 PetscErrorCode ierr; 2141 2142 PetscFunctionBegin; 2143 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2144 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2145 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2146 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2147 PetscFunctionReturn(0); 2148 } 2149 2150 #undef __FUNCT__ 2151 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2152 /*@C 2153 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2154 2155 Collective on MPI_Comm 2156 2157 Input Parameters: 2158 + A - the matrix 2159 - data - the array (or NULL) 2160 2161 Notes: 2162 The data input variable is intended primarily for Fortran programmers 2163 who wish to allocate their own matrix memory space. Most users should 2164 need not call this routine. 2165 2166 Level: intermediate 2167 2168 .keywords: dense, matrix, LAPACK, BLAS 2169 2170 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2171 2172 @*/ 2173 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2174 { 2175 PetscErrorCode ierr; 2176 2177 PetscFunctionBegin; 2178 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2179 PetscFunctionReturn(0); 2180 } 2181 2182 #undef __FUNCT__ 2183 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2184 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2185 { 2186 Mat_SeqDense *b; 2187 PetscErrorCode ierr; 2188 2189 PetscFunctionBegin; 2190 B->preallocated = PETSC_TRUE; 2191 2192 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2193 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2194 2195 b = (Mat_SeqDense*)B->data; 2196 b->Mmax = B->rmap->n; 2197 b->Nmax = B->cmap->n; 2198 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2199 2200 if (!data) { /* petsc-allocated storage */ 2201 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2202 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2203 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2204 2205 b->user_alloc = PETSC_FALSE; 2206 } else { /* user-allocated storage */ 2207 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2208 b->v = data; 2209 b->user_alloc = PETSC_TRUE; 2210 } 2211 B->assembled = PETSC_TRUE; 2212 PetscFunctionReturn(0); 2213 } 2214 2215 #undef __FUNCT__ 2216 #define __FUNCT__ "MatSeqDenseSetLDA" 2217 /*@C 2218 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2219 2220 Input parameter: 2221 + A - the matrix 2222 - lda - the leading dimension 2223 2224 Notes: 2225 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2226 it asserts that the preallocation has a leading dimension (the LDA parameter 2227 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2228 2229 Level: intermediate 2230 2231 .keywords: dense, matrix, LAPACK, BLAS 2232 2233 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2234 2235 @*/ 2236 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2237 { 2238 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2239 2240 PetscFunctionBegin; 2241 if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2242 b->lda = lda; 2243 b->changelda = PETSC_FALSE; 2244 b->Mmax = PetscMax(b->Mmax,lda); 2245 PetscFunctionReturn(0); 2246 } 2247 2248 /*MC 2249 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2250 2251 Options Database Keys: 2252 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2253 2254 Level: beginner 2255 2256 .seealso: MatCreateSeqDense() 2257 2258 M*/ 2259 2260 #undef __FUNCT__ 2261 #define __FUNCT__ "MatCreate_SeqDense" 2262 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2263 { 2264 Mat_SeqDense *b; 2265 PetscErrorCode ierr; 2266 PetscMPIInt size; 2267 2268 PetscFunctionBegin; 2269 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2270 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2271 2272 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2273 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2274 B->data = (void*)b; 2275 2276 b->pivots = 0; 2277 b->roworiented = PETSC_TRUE; 2278 b->v = 0; 2279 b->changelda = PETSC_FALSE; 2280 2281 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2282 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2283 2284 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2285 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2288 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2289 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2292 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2293 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2294 PetscFunctionReturn(0); 2295 } 2296