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 = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatGetInfo_SeqDense" 79 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 80 { 81 PetscInt N = A->rmap->n*A->cmap->n; 82 83 PetscFunctionBegin; 84 info->block_size = 1.0; 85 info->nz_allocated = (double)N; 86 info->nz_used = (double)N; 87 info->nz_unneeded = (double)0; 88 info->assemblies = (double)A->num_ass; 89 info->mallocs = 0; 90 info->memory = ((PetscObject)A)->mem; 91 info->fill_ratio_given = 0; 92 info->fill_ratio_needed = 0; 93 info->factor_mallocs = 0; 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatScale_SeqDense" 99 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 100 { 101 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 102 PetscScalar oalpha = alpha; 103 PetscErrorCode ierr; 104 PetscBLASInt one = 1,j,nz,lda; 105 106 PetscFunctionBegin; 107 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 108 if (lda>A->rmap->n) { 109 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 110 for (j=0; j<A->cmap->n; j++) { 111 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one)); 112 } 113 } else { 114 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 115 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one)); 116 } 117 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatIsHermitian_SeqDense" 123 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 124 { 125 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 126 PetscInt i,j,m = A->rmap->n,N; 127 PetscScalar *v = a->v; 128 129 PetscFunctionBegin; 130 *fl = PETSC_FALSE; 131 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 132 N = a->lda; 133 134 for (i=0; i<m; i++) { 135 for (j=i+1; j<m; j++) { 136 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 137 } 138 } 139 *fl = PETSC_TRUE; 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 145 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 146 { 147 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 148 PetscErrorCode ierr; 149 PetscInt lda = (PetscInt)mat->lda,j,m; 150 151 PetscFunctionBegin; 152 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 153 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 154 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 155 if (cpvalues == MAT_COPY_VALUES) { 156 l = (Mat_SeqDense*)newi->data; 157 if (lda>A->rmap->n) { 158 m = A->rmap->n; 159 for (j=0; j<A->cmap->n; j++) { 160 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 161 } 162 } else { 163 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 164 } 165 } 166 newi->assembled = PETSC_TRUE; 167 PetscFunctionReturn(0); 168 } 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "MatDuplicate_SeqDense" 172 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 173 { 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 178 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 179 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 180 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 185 extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 189 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 190 { 191 MatFactorInfo info; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 196 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatSolve_SeqDense" 202 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 203 { 204 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 205 PetscErrorCode ierr; 206 PetscScalar *x,*y; 207 PetscBLASInt one = 1,info,m; 208 209 PetscFunctionBegin; 210 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 211 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 212 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 213 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 214 if (A->factortype == MAT_FACTOR_LU) { 215 #if defined(PETSC_MISSING_LAPACK_GETRS) 216 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 217 #else 218 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 219 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 220 #endif 221 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 222 #if defined(PETSC_MISSING_LAPACK_POTRS) 223 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 224 #else 225 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 226 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 227 #endif 228 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 229 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 230 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 231 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "MatMatSolve_SeqDense" 237 PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 238 { 239 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 240 PetscErrorCode ierr; 241 PetscScalar *b,*x; 242 PetscInt n; 243 PetscBLASInt nrhs,info,m; 244 PetscBool flg; 245 246 PetscFunctionBegin; 247 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 248 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 249 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 250 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 251 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 252 253 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 254 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 255 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 256 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 257 258 ierr = PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));CHKERRQ(ierr); 259 260 if (A->factortype == MAT_FACTOR_LU) { 261 #if defined(PETSC_MISSING_LAPACK_GETRS) 262 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 263 #else 264 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 265 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 266 #endif 267 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 268 #if defined(PETSC_MISSING_LAPACK_POTRS) 269 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 270 #else 271 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 272 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 273 #endif 274 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 275 276 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 277 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 278 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatSolveTranspose_SeqDense" 284 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 285 { 286 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 287 PetscErrorCode ierr; 288 PetscScalar *x,*y; 289 PetscBLASInt one = 1,info,m; 290 291 PetscFunctionBegin; 292 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 293 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 294 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 295 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 296 /* assume if pivots exist then use LU; else Cholesky */ 297 if (mat->pivots) { 298 #if defined(PETSC_MISSING_LAPACK_GETRS) 299 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 300 #else 301 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 302 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 303 #endif 304 } else { 305 #if defined(PETSC_MISSING_LAPACK_POTRS) 306 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 307 #else 308 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 309 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 310 #endif 311 } 312 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 313 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 314 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 315 PetscFunctionReturn(0); 316 } 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatSolveAdd_SeqDense" 320 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 321 { 322 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 323 PetscErrorCode ierr; 324 PetscScalar *x,*y,sone = 1.0; 325 Vec tmp = 0; 326 PetscBLASInt one = 1,info,m; 327 328 PetscFunctionBegin; 329 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 330 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 331 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 332 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 333 if (yy == zz) { 334 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 335 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 336 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 337 } 338 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 339 /* assume if pivots exist then use LU; else Cholesky */ 340 if (mat->pivots) { 341 #if defined(PETSC_MISSING_LAPACK_GETRS) 342 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 343 #else 344 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 345 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 346 #endif 347 } else { 348 #if defined(PETSC_MISSING_LAPACK_POTRS) 349 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 350 #else 351 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 352 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 353 #endif 354 } 355 if (tmp) { 356 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 357 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 358 } else { 359 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 360 } 361 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 362 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 363 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 364 PetscFunctionReturn(0); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 369 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 370 { 371 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 372 PetscErrorCode ierr; 373 PetscScalar *x,*y,sone = 1.0; 374 Vec tmp; 375 PetscBLASInt one = 1,info,m; 376 377 PetscFunctionBegin; 378 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 379 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 380 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 381 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 382 if (yy == zz) { 383 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 384 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 385 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 386 } 387 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 388 /* assume if pivots exist then use LU; else Cholesky */ 389 if (mat->pivots) { 390 #if defined(PETSC_MISSING_LAPACK_GETRS) 391 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 392 #else 393 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 394 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 395 #endif 396 } else { 397 #if defined(PETSC_MISSING_LAPACK_POTRS) 398 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 399 #else 400 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 401 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 402 #endif 403 } 404 if (tmp) { 405 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 406 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 407 } else { 408 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 409 } 410 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 411 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 412 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 413 PetscFunctionReturn(0); 414 } 415 416 /* ---------------------------------------------------------------*/ 417 /* COMMENT: I have chosen to hide row permutation in the pivots, 418 rather than put it in the Mat->row slot.*/ 419 #undef __FUNCT__ 420 #define __FUNCT__ "MatLUFactor_SeqDense" 421 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 422 { 423 #if defined(PETSC_MISSING_LAPACK_GETRF) 424 PetscFunctionBegin; 425 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 426 #else 427 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 428 PetscErrorCode ierr; 429 PetscBLASInt n,m,info; 430 431 PetscFunctionBegin; 432 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 433 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 434 if (!mat->pivots) { 435 ierr = PetscMalloc1((A->rmap->n+1),&mat->pivots);CHKERRQ(ierr); 436 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 437 } 438 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 439 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 440 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 441 ierr = PetscFPTrapPop();CHKERRQ(ierr); 442 443 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 444 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 445 A->ops->solve = MatSolve_SeqDense; 446 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 447 A->ops->solveadd = MatSolveAdd_SeqDense; 448 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 449 A->factortype = MAT_FACTOR_LU; 450 451 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 452 #endif 453 PetscFunctionReturn(0); 454 } 455 456 #undef __FUNCT__ 457 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 458 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 459 { 460 #if defined(PETSC_MISSING_LAPACK_POTRF) 461 PetscFunctionBegin; 462 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 463 #else 464 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 465 PetscErrorCode ierr; 466 PetscBLASInt info,n; 467 468 PetscFunctionBegin; 469 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 470 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 471 472 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 473 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 474 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 475 A->ops->solve = MatSolve_SeqDense; 476 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 477 A->ops->solveadd = MatSolveAdd_SeqDense; 478 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 479 A->factortype = MAT_FACTOR_CHOLESKY; 480 481 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 482 #endif 483 PetscFunctionReturn(0); 484 } 485 486 487 #undef __FUNCT__ 488 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 489 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 490 { 491 PetscErrorCode ierr; 492 MatFactorInfo info; 493 494 PetscFunctionBegin; 495 info.fill = 1.0; 496 497 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 498 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 499 PetscFunctionReturn(0); 500 } 501 502 #undef __FUNCT__ 503 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 504 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 505 { 506 PetscFunctionBegin; 507 fact->assembled = PETSC_TRUE; 508 fact->preallocated = PETSC_TRUE; 509 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 510 PetscFunctionReturn(0); 511 } 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 515 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 516 { 517 PetscFunctionBegin; 518 fact->preallocated = PETSC_TRUE; 519 fact->assembled = PETSC_TRUE; 520 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 521 PetscFunctionReturn(0); 522 } 523 524 #undef __FUNCT__ 525 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 526 PETSC_EXTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 527 { 528 PetscErrorCode ierr; 529 530 PetscFunctionBegin; 531 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 532 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 533 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 534 if (ftype == MAT_FACTOR_LU) { 535 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 536 } else { 537 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 538 } 539 (*fact)->factortype = ftype; 540 PetscFunctionReturn(0); 541 } 542 543 /* ------------------------------------------------------------------*/ 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatSOR_SeqDense" 546 PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 547 { 548 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 549 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 550 PetscErrorCode ierr; 551 PetscInt m = A->rmap->n,i; 552 PetscBLASInt o = 1,bm; 553 554 PetscFunctionBegin; 555 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 556 if (flag & SOR_ZERO_INITIAL_GUESS) { 557 /* this is a hack fix, should have another version without the second BLASdot */ 558 ierr = VecSet(xx,zero);CHKERRQ(ierr); 559 } 560 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 561 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 562 its = its*lits; 563 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 564 while (its--) { 565 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 566 for (i=0; i<m; i++) { 567 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 568 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 569 } 570 } 571 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 572 for (i=m-1; i>=0; i--) { 573 PetscStackCallBLAS("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 574 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 575 } 576 } 577 } 578 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 579 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 /* -----------------------------------------------------------------*/ 584 #undef __FUNCT__ 585 #define __FUNCT__ "MatMultTranspose_SeqDense" 586 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 587 { 588 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 589 PetscScalar *v = mat->v,*x,*y; 590 PetscErrorCode ierr; 591 PetscBLASInt m, n,_One=1; 592 PetscScalar _DOne=1.0,_DZero=0.0; 593 594 PetscFunctionBegin; 595 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 596 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 597 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 598 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 599 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 600 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 601 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 602 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 603 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607 #undef __FUNCT__ 608 #define __FUNCT__ "MatMult_SeqDense" 609 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 610 { 611 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 612 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 613 PetscErrorCode ierr; 614 PetscBLASInt m, n, _One=1; 615 616 PetscFunctionBegin; 617 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 618 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 619 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 620 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 621 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 622 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 623 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 624 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 625 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "MatMultAdd_SeqDense" 631 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 632 { 633 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 634 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 635 PetscErrorCode ierr; 636 PetscBLASInt m, n, _One=1; 637 638 PetscFunctionBegin; 639 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 640 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 641 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 642 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 643 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 644 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 645 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 646 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 647 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 648 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 #undef __FUNCT__ 653 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 654 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 655 { 656 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 657 PetscScalar *v = mat->v,*x,*y; 658 PetscErrorCode ierr; 659 PetscBLASInt m, n, _One=1; 660 PetscScalar _DOne=1.0; 661 662 PetscFunctionBegin; 663 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 664 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 665 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 666 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 667 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 668 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 669 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 670 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 671 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 672 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 673 PetscFunctionReturn(0); 674 } 675 676 /* -----------------------------------------------------------------*/ 677 #undef __FUNCT__ 678 #define __FUNCT__ "MatGetRow_SeqDense" 679 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 680 { 681 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 682 PetscScalar *v; 683 PetscErrorCode ierr; 684 PetscInt i; 685 686 PetscFunctionBegin; 687 *ncols = A->cmap->n; 688 if (cols) { 689 ierr = PetscMalloc1((A->cmap->n+1),cols);CHKERRQ(ierr); 690 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 691 } 692 if (vals) { 693 ierr = PetscMalloc1((A->cmap->n+1),vals);CHKERRQ(ierr); 694 v = mat->v + row; 695 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 696 } 697 PetscFunctionReturn(0); 698 } 699 700 #undef __FUNCT__ 701 #define __FUNCT__ "MatRestoreRow_SeqDense" 702 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 703 { 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 708 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 709 PetscFunctionReturn(0); 710 } 711 /* ----------------------------------------------------------------*/ 712 #undef __FUNCT__ 713 #define __FUNCT__ "MatSetValues_SeqDense" 714 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 715 { 716 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 717 PetscInt i,j,idx=0; 718 719 PetscFunctionBegin; 720 if (v) PetscValidScalarPointer(v,6); 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