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 (!mat->roworiented) { 721 if (addv == INSERT_VALUES) { 722 for (j=0; j<n; j++) { 723 if (indexn[j] < 0) {idx += m; continue;} 724 #if defined(PETSC_USE_DEBUG) 725 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); 726 #endif 727 for (i=0; i<m; i++) { 728 if (indexm[i] < 0) {idx++; continue;} 729 #if defined(PETSC_USE_DEBUG) 730 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); 731 #endif 732 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 733 } 734 } 735 } else { 736 for (j=0; j<n; j++) { 737 if (indexn[j] < 0) {idx += m; continue;} 738 #if defined(PETSC_USE_DEBUG) 739 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); 740 #endif 741 for (i=0; i<m; i++) { 742 if (indexm[i] < 0) {idx++; continue;} 743 #if defined(PETSC_USE_DEBUG) 744 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); 745 #endif 746 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 747 } 748 } 749 } 750 } else { 751 if (addv == INSERT_VALUES) { 752 for (i=0; i<m; i++) { 753 if (indexm[i] < 0) { idx += n; continue;} 754 #if defined(PETSC_USE_DEBUG) 755 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); 756 #endif 757 for (j=0; j<n; j++) { 758 if (indexn[j] < 0) { idx++; continue;} 759 #if defined(PETSC_USE_DEBUG) 760 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); 761 #endif 762 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 763 } 764 } 765 } else { 766 for (i=0; i<m; i++) { 767 if (indexm[i] < 0) { idx += n; continue;} 768 #if defined(PETSC_USE_DEBUG) 769 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); 770 #endif 771 for (j=0; j<n; j++) { 772 if (indexn[j] < 0) { idx++; continue;} 773 #if defined(PETSC_USE_DEBUG) 774 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); 775 #endif 776 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 777 } 778 } 779 } 780 } 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "MatGetValues_SeqDense" 786 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 787 { 788 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 789 PetscInt i,j; 790 791 PetscFunctionBegin; 792 /* row-oriented output */ 793 for (i=0; i<m; i++) { 794 if (indexm[i] < 0) {v += n;continue;} 795 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); 796 for (j=0; j<n; j++) { 797 if (indexn[j] < 0) {v++; continue;} 798 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); 799 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 800 } 801 } 802 PetscFunctionReturn(0); 803 } 804 805 /* -----------------------------------------------------------------*/ 806 807 #undef __FUNCT__ 808 #define __FUNCT__ "MatLoad_SeqDense" 809 PetscErrorCode MatLoad_SeqDense(Mat newmat,PetscViewer viewer) 810 { 811 Mat_SeqDense *a; 812 PetscErrorCode ierr; 813 PetscInt *scols,i,j,nz,header[4]; 814 int fd; 815 PetscMPIInt size; 816 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 817 PetscScalar *vals,*svals,*v,*w; 818 MPI_Comm comm; 819 820 PetscFunctionBegin; 821 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 822 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 823 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 824 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 825 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 826 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 827 M = header[1]; N = header[2]; nz = header[3]; 828 829 /* set global size if not set already*/ 830 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 831 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 832 } else { 833 /* if sizes and type are already set, check if the vector global sizes are correct */ 834 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 835 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); 836 } 837 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 838 839 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 840 a = (Mat_SeqDense*)newmat->data; 841 v = a->v; 842 /* Allocate some temp space to read in the values and then flip them 843 from row major to column major */ 844 ierr = PetscMalloc1((M*N > 0 ? M*N : 1),&w);CHKERRQ(ierr); 845 /* read in nonzero values */ 846 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 847 /* now flip the values and store them in the matrix*/ 848 for (j=0; j<N; j++) { 849 for (i=0; i<M; i++) { 850 *v++ =w[i*N+j]; 851 } 852 } 853 ierr = PetscFree(w);CHKERRQ(ierr); 854 } else { 855 /* read row lengths */ 856 ierr = PetscMalloc1((M+1),&rowlengths);CHKERRQ(ierr); 857 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 858 859 a = (Mat_SeqDense*)newmat->data; 860 v = a->v; 861 862 /* read column indices and nonzeros */ 863 ierr = PetscMalloc1((nz+1),&scols);CHKERRQ(ierr); 864 cols = scols; 865 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 866 ierr = PetscMalloc1((nz+1),&svals);CHKERRQ(ierr); 867 vals = svals; 868 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 869 870 /* insert into matrix */ 871 for (i=0; i<M; i++) { 872 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 873 svals += rowlengths[i]; scols += rowlengths[i]; 874 } 875 ierr = PetscFree(vals);CHKERRQ(ierr); 876 ierr = PetscFree(cols);CHKERRQ(ierr); 877 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 878 } 879 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 880 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 881 PetscFunctionReturn(0); 882 } 883 884 #undef __FUNCT__ 885 #define __FUNCT__ "MatView_SeqDense_ASCII" 886 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 887 { 888 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 889 PetscErrorCode ierr; 890 PetscInt i,j; 891 const char *name; 892 PetscScalar *v; 893 PetscViewerFormat format; 894 #if defined(PETSC_USE_COMPLEX) 895 PetscBool allreal = PETSC_TRUE; 896 #endif 897 898 PetscFunctionBegin; 899 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 900 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 901 PetscFunctionReturn(0); /* do nothing for now */ 902 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 903 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 904 for (i=0; i<A->rmap->n; i++) { 905 v = a->v + i; 906 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 907 for (j=0; j<A->cmap->n; j++) { 908 #if defined(PETSC_USE_COMPLEX) 909 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 910 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 911 } else if (PetscRealPart(*v)) { 912 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 913 } 914 #else 915 if (*v) { 916 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 917 } 918 #endif 919 v += a->lda; 920 } 921 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 922 } 923 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 924 } else { 925 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 926 #if defined(PETSC_USE_COMPLEX) 927 /* determine if matrix has all real values */ 928 v = a->v; 929 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 930 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 931 } 932 #endif 933 if (format == PETSC_VIEWER_ASCII_MATLAB) { 934 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 935 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 936 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 937 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 938 } 939 940 for (i=0; i<A->rmap->n; i++) { 941 v = a->v + i; 942 for (j=0; j<A->cmap->n; j++) { 943 #if defined(PETSC_USE_COMPLEX) 944 if (allreal) { 945 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 946 } else { 947 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 948 } 949 #else 950 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 951 #endif 952 v += a->lda; 953 } 954 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 955 } 956 if (format == PETSC_VIEWER_ASCII_MATLAB) { 957 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 958 } 959 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 960 } 961 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "MatView_SeqDense_Binary" 967 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 968 { 969 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 970 PetscErrorCode ierr; 971 int fd; 972 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 973 PetscScalar *v,*anonz,*vals; 974 PetscViewerFormat format; 975 976 PetscFunctionBegin; 977 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 978 979 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 980 if (format == PETSC_VIEWER_NATIVE) { 981 /* store the matrix as a dense matrix */ 982 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 983 984 col_lens[0] = MAT_FILE_CLASSID; 985 col_lens[1] = m; 986 col_lens[2] = n; 987 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 988 989 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 990 ierr = PetscFree(col_lens);CHKERRQ(ierr); 991 992 /* write out matrix, by rows */ 993 ierr = PetscMalloc1((m*n+1),&vals);CHKERRQ(ierr); 994 v = a->v; 995 for (j=0; j<n; j++) { 996 for (i=0; i<m; i++) { 997 vals[j + i*n] = *v++; 998 } 999 } 1000 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1001 ierr = PetscFree(vals);CHKERRQ(ierr); 1002 } else { 1003 ierr = PetscMalloc1((4+nz),&col_lens);CHKERRQ(ierr); 1004 1005 col_lens[0] = MAT_FILE_CLASSID; 1006 col_lens[1] = m; 1007 col_lens[2] = n; 1008 col_lens[3] = nz; 1009 1010 /* store lengths of each row and write (including header) to file */ 1011 for (i=0; i<m; i++) col_lens[4+i] = n; 1012 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1013 1014 /* Possibly should write in smaller increments, not whole matrix at once? */ 1015 /* store column indices (zero start index) */ 1016 ict = 0; 1017 for (i=0; i<m; i++) { 1018 for (j=0; j<n; j++) col_lens[ict++] = j; 1019 } 1020 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1021 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1022 1023 /* store nonzero values */ 1024 ierr = PetscMalloc1((nz+1),&anonz);CHKERRQ(ierr); 1025 ict = 0; 1026 for (i=0; i<m; i++) { 1027 v = a->v + i; 1028 for (j=0; j<n; j++) { 1029 anonz[ict++] = *v; v += a->lda; 1030 } 1031 } 1032 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1033 ierr = PetscFree(anonz);CHKERRQ(ierr); 1034 } 1035 PetscFunctionReturn(0); 1036 } 1037 1038 #include <petscdraw.h> 1039 #undef __FUNCT__ 1040 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1041 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1042 { 1043 Mat A = (Mat) Aa; 1044 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1045 PetscErrorCode ierr; 1046 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1047 PetscScalar *v = a->v; 1048 PetscViewer viewer; 1049 PetscDraw popup; 1050 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1051 PetscViewerFormat format; 1052 1053 PetscFunctionBegin; 1054 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1055 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1056 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1057 1058 /* Loop over matrix elements drawing boxes */ 1059 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1060 /* Blue for negative and Red for positive */ 1061 color = PETSC_DRAW_BLUE; 1062 for (j = 0; j < n; j++) { 1063 x_l = j; 1064 x_r = x_l + 1.0; 1065 for (i = 0; i < m; i++) { 1066 y_l = m - i - 1.0; 1067 y_r = y_l + 1.0; 1068 if (PetscRealPart(v[j*m+i]) > 0.) { 1069 color = PETSC_DRAW_RED; 1070 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1071 color = PETSC_DRAW_BLUE; 1072 } else { 1073 continue; 1074 } 1075 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1076 } 1077 } 1078 } else { 1079 /* use contour shading to indicate magnitude of values */ 1080 /* first determine max of all nonzero values */ 1081 for (i = 0; i < m*n; i++) { 1082 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1083 } 1084 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1085 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1086 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1087 for (j = 0; j < n; j++) { 1088 x_l = j; 1089 x_r = x_l + 1.0; 1090 for (i = 0; i < m; i++) { 1091 y_l = m - i - 1.0; 1092 y_r = y_l + 1.0; 1093 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1094 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1095 } 1096 } 1097 } 1098 PetscFunctionReturn(0); 1099 } 1100 1101 #undef __FUNCT__ 1102 #define __FUNCT__ "MatView_SeqDense_Draw" 1103 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1104 { 1105 PetscDraw draw; 1106 PetscBool isnull; 1107 PetscReal xr,yr,xl,yl,h,w; 1108 PetscErrorCode ierr; 1109 1110 PetscFunctionBegin; 1111 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1112 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1113 if (isnull) PetscFunctionReturn(0); 1114 1115 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1116 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1117 xr += w; yr += h; xl = -w; yl = -h; 1118 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1119 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1120 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 #undef __FUNCT__ 1125 #define __FUNCT__ "MatView_SeqDense" 1126 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1127 { 1128 PetscErrorCode ierr; 1129 PetscBool iascii,isbinary,isdraw; 1130 1131 PetscFunctionBegin; 1132 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1134 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1135 1136 if (iascii) { 1137 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1138 } else if (isbinary) { 1139 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1140 } else if (isdraw) { 1141 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1142 } 1143 PetscFunctionReturn(0); 1144 } 1145 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "MatDestroy_SeqDense" 1148 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1149 { 1150 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1151 PetscErrorCode ierr; 1152 1153 PetscFunctionBegin; 1154 #if defined(PETSC_USE_LOG) 1155 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1156 #endif 1157 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1158 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1159 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1160 1161 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1162 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1163 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1164 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatTranspose_SeqDense" 1176 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1177 { 1178 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1179 PetscErrorCode ierr; 1180 PetscInt k,j,m,n,M; 1181 PetscScalar *v,tmp; 1182 1183 PetscFunctionBegin; 1184 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1185 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1186 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1187 else { 1188 for (j=0; j<m; j++) { 1189 for (k=0; k<j; k++) { 1190 tmp = v[j + k*M]; 1191 v[j + k*M] = v[k + j*M]; 1192 v[k + j*M] = tmp; 1193 } 1194 } 1195 } 1196 } else { /* out-of-place transpose */ 1197 Mat tmat; 1198 Mat_SeqDense *tmatd; 1199 PetscScalar *v2; 1200 1201 if (reuse == MAT_INITIAL_MATRIX) { 1202 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1203 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1204 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1205 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1206 } else { 1207 tmat = *matout; 1208 } 1209 tmatd = (Mat_SeqDense*)tmat->data; 1210 v = mat->v; v2 = tmatd->v; 1211 for (j=0; j<n; j++) { 1212 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1213 } 1214 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1215 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1216 *matout = tmat; 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatEqual_SeqDense" 1223 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1224 { 1225 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1226 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1227 PetscInt i,j; 1228 PetscScalar *v1,*v2; 1229 1230 PetscFunctionBegin; 1231 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1232 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1233 for (i=0; i<A1->rmap->n; i++) { 1234 v1 = mat1->v+i; v2 = mat2->v+i; 1235 for (j=0; j<A1->cmap->n; j++) { 1236 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1237 v1 += mat1->lda; v2 += mat2->lda; 1238 } 1239 } 1240 *flg = PETSC_TRUE; 1241 PetscFunctionReturn(0); 1242 } 1243 1244 #undef __FUNCT__ 1245 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1246 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1247 { 1248 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1249 PetscErrorCode ierr; 1250 PetscInt i,n,len; 1251 PetscScalar *x,zero = 0.0; 1252 1253 PetscFunctionBegin; 1254 ierr = VecSet(v,zero);CHKERRQ(ierr); 1255 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1256 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1257 len = PetscMin(A->rmap->n,A->cmap->n); 1258 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1259 for (i=0; i<len; i++) { 1260 x[i] = mat->v[i*mat->lda + i]; 1261 } 1262 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 #undef __FUNCT__ 1267 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1268 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1269 { 1270 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1271 PetscScalar *l,*r,x,*v; 1272 PetscErrorCode ierr; 1273 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1274 1275 PetscFunctionBegin; 1276 if (ll) { 1277 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1278 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1279 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1280 for (i=0; i<m; i++) { 1281 x = l[i]; 1282 v = mat->v + i; 1283 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1284 } 1285 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1286 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1287 } 1288 if (rr) { 1289 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1290 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1291 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1292 for (i=0; i<n; i++) { 1293 x = r[i]; 1294 v = mat->v + i*m; 1295 for (j=0; j<m; j++) (*v++) *= x; 1296 } 1297 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1298 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1299 } 1300 PetscFunctionReturn(0); 1301 } 1302 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatNorm_SeqDense" 1305 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1306 { 1307 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1308 PetscScalar *v = mat->v; 1309 PetscReal sum = 0.0; 1310 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1311 PetscErrorCode ierr; 1312 1313 PetscFunctionBegin; 1314 if (type == NORM_FROBENIUS) { 1315 if (lda>m) { 1316 for (j=0; j<A->cmap->n; j++) { 1317 v = mat->v+j*lda; 1318 for (i=0; i<m; i++) { 1319 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1320 } 1321 } 1322 } else { 1323 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1324 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1325 } 1326 } 1327 *nrm = PetscSqrtReal(sum); 1328 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1329 } else if (type == NORM_1) { 1330 *nrm = 0.0; 1331 for (j=0; j<A->cmap->n; j++) { 1332 v = mat->v + j*mat->lda; 1333 sum = 0.0; 1334 for (i=0; i<A->rmap->n; i++) { 1335 sum += PetscAbsScalar(*v); v++; 1336 } 1337 if (sum > *nrm) *nrm = sum; 1338 } 1339 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1340 } else if (type == NORM_INFINITY) { 1341 *nrm = 0.0; 1342 for (j=0; j<A->rmap->n; j++) { 1343 v = mat->v + j; 1344 sum = 0.0; 1345 for (i=0; i<A->cmap->n; i++) { 1346 sum += PetscAbsScalar(*v); v += mat->lda; 1347 } 1348 if (sum > *nrm) *nrm = sum; 1349 } 1350 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1351 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1352 PetscFunctionReturn(0); 1353 } 1354 1355 #undef __FUNCT__ 1356 #define __FUNCT__ "MatSetOption_SeqDense" 1357 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1358 { 1359 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1360 PetscErrorCode ierr; 1361 1362 PetscFunctionBegin; 1363 switch (op) { 1364 case MAT_ROW_ORIENTED: 1365 aij->roworiented = flg; 1366 break; 1367 case MAT_NEW_NONZERO_LOCATIONS: 1368 case MAT_NEW_NONZERO_LOCATION_ERR: 1369 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1370 case MAT_NEW_DIAGONALS: 1371 case MAT_KEEP_NONZERO_PATTERN: 1372 case MAT_IGNORE_OFF_PROC_ENTRIES: 1373 case MAT_USE_HASH_TABLE: 1374 case MAT_IGNORE_LOWER_TRIANGULAR: 1375 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1376 break; 1377 case MAT_SPD: 1378 case MAT_SYMMETRIC: 1379 case MAT_STRUCTURALLY_SYMMETRIC: 1380 case MAT_HERMITIAN: 1381 case MAT_SYMMETRY_ETERNAL: 1382 /* These options are handled directly by MatSetOption() */ 1383 break; 1384 default: 1385 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1386 } 1387 PetscFunctionReturn(0); 1388 } 1389 1390 #undef __FUNCT__ 1391 #define __FUNCT__ "MatZeroEntries_SeqDense" 1392 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1393 { 1394 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1395 PetscErrorCode ierr; 1396 PetscInt lda=l->lda,m=A->rmap->n,j; 1397 1398 PetscFunctionBegin; 1399 if (lda>m) { 1400 for (j=0; j<A->cmap->n; j++) { 1401 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1402 } 1403 } else { 1404 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1405 } 1406 PetscFunctionReturn(0); 1407 } 1408 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "MatZeroRows_SeqDense" 1411 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1412 { 1413 PetscErrorCode ierr; 1414 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1415 PetscInt m = l->lda, n = A->cmap->n, i,j; 1416 PetscScalar *slot,*bb; 1417 const PetscScalar *xx; 1418 1419 PetscFunctionBegin; 1420 #if defined(PETSC_USE_DEBUG) 1421 for (i=0; i<N; i++) { 1422 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1423 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); 1424 } 1425 #endif 1426 1427 /* fix right hand side if needed */ 1428 if (x && b) { 1429 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1430 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1431 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1432 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1433 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1434 } 1435 1436 for (i=0; i<N; i++) { 1437 slot = l->v + rows[i]; 1438 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1439 } 1440 if (diag != 0.0) { 1441 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1442 for (i=0; i<N; i++) { 1443 slot = l->v + (m+1)*rows[i]; 1444 *slot = diag; 1445 } 1446 } 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "MatDenseGetArray_SeqDense" 1452 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1453 { 1454 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1455 1456 PetscFunctionBegin; 1457 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"); 1458 *array = mat->v; 1459 PetscFunctionReturn(0); 1460 } 1461 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatDenseRestoreArray_SeqDense" 1464 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1465 { 1466 PetscFunctionBegin; 1467 *array = 0; /* user cannot accidently use the array later */ 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "MatDenseGetArray" 1473 /*@C 1474 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1475 1476 Not Collective 1477 1478 Input Parameter: 1479 . mat - a MATSEQDENSE matrix 1480 1481 Output Parameter: 1482 . array - pointer to the data 1483 1484 Level: intermediate 1485 1486 .seealso: MatDenseRestoreArray() 1487 @*/ 1488 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1489 { 1490 PetscErrorCode ierr; 1491 1492 PetscFunctionBegin; 1493 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1494 PetscFunctionReturn(0); 1495 } 1496 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "MatDenseRestoreArray" 1499 /*@C 1500 MatDenseRestoreArray - returns access to the array where the data for a SeqDense matrix is stored obtained by MatDenseGetArray() 1501 1502 Not Collective 1503 1504 Input Parameters: 1505 . mat - a MATSEQDENSE matrix 1506 . array - pointer to the data 1507 1508 Level: intermediate 1509 1510 .seealso: MatDenseGetArray() 1511 @*/ 1512 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1513 { 1514 PetscErrorCode ierr; 1515 1516 PetscFunctionBegin; 1517 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1518 PetscFunctionReturn(0); 1519 } 1520 1521 #undef __FUNCT__ 1522 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1523 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1524 { 1525 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1526 PetscErrorCode ierr; 1527 PetscInt i,j,nrows,ncols; 1528 const PetscInt *irow,*icol; 1529 PetscScalar *av,*bv,*v = mat->v; 1530 Mat newmat; 1531 1532 PetscFunctionBegin; 1533 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1534 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1535 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1536 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1537 1538 /* Check submatrixcall */ 1539 if (scall == MAT_REUSE_MATRIX) { 1540 PetscInt n_cols,n_rows; 1541 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1542 if (n_rows != nrows || n_cols != ncols) { 1543 /* resize the result matrix to match number of requested rows/columns */ 1544 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1545 } 1546 newmat = *B; 1547 } else { 1548 /* Create and fill new matrix */ 1549 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1550 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1551 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1552 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1553 } 1554 1555 /* Now extract the data pointers and do the copy,column at a time */ 1556 bv = ((Mat_SeqDense*)newmat->data)->v; 1557 1558 for (i=0; i<ncols; i++) { 1559 av = v + mat->lda*icol[i]; 1560 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1561 } 1562 1563 /* Assemble the matrices so that the correct flags are set */ 1564 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1565 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1566 1567 /* Free work space */ 1568 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1569 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1570 *B = newmat; 1571 PetscFunctionReturn(0); 1572 } 1573 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1576 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1577 { 1578 PetscErrorCode ierr; 1579 PetscInt i; 1580 1581 PetscFunctionBegin; 1582 if (scall == MAT_INITIAL_MATRIX) { 1583 ierr = PetscMalloc1((n+1),B);CHKERRQ(ierr); 1584 } 1585 1586 for (i=0; i<n; i++) { 1587 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1588 } 1589 PetscFunctionReturn(0); 1590 } 1591 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1594 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1595 { 1596 PetscFunctionBegin; 1597 PetscFunctionReturn(0); 1598 } 1599 1600 #undef __FUNCT__ 1601 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1602 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1603 { 1604 PetscFunctionBegin; 1605 PetscFunctionReturn(0); 1606 } 1607 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "MatCopy_SeqDense" 1610 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1611 { 1612 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 1613 PetscErrorCode ierr; 1614 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1615 1616 PetscFunctionBegin; 1617 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1618 if (A->ops->copy != B->ops->copy) { 1619 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1620 PetscFunctionReturn(0); 1621 } 1622 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1623 if (lda1>m || lda2>m) { 1624 for (j=0; j<n; j++) { 1625 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1626 } 1627 } else { 1628 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1629 } 1630 PetscFunctionReturn(0); 1631 } 1632 1633 #undef __FUNCT__ 1634 #define __FUNCT__ "MatSetUp_SeqDense" 1635 PetscErrorCode MatSetUp_SeqDense(Mat A) 1636 { 1637 PetscErrorCode ierr; 1638 1639 PetscFunctionBegin; 1640 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "MatConjugate_SeqDense" 1646 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1647 { 1648 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1649 PetscInt i,nz = A->rmap->n*A->cmap->n; 1650 PetscScalar *aa = a->v; 1651 1652 PetscFunctionBegin; 1653 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "MatRealPart_SeqDense" 1659 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1660 { 1661 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1662 PetscInt i,nz = A->rmap->n*A->cmap->n; 1663 PetscScalar *aa = a->v; 1664 1665 PetscFunctionBegin; 1666 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1667 PetscFunctionReturn(0); 1668 } 1669 1670 #undef __FUNCT__ 1671 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1672 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1673 { 1674 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1675 PetscInt i,nz = A->rmap->n*A->cmap->n; 1676 PetscScalar *aa = a->v; 1677 1678 PetscFunctionBegin; 1679 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1680 PetscFunctionReturn(0); 1681 } 1682 1683 /* ----------------------------------------------------------------*/ 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1686 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1687 { 1688 PetscErrorCode ierr; 1689 1690 PetscFunctionBegin; 1691 if (scall == MAT_INITIAL_MATRIX) { 1692 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1693 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1694 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1695 } 1696 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1697 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1698 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1704 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1705 { 1706 PetscErrorCode ierr; 1707 PetscInt m=A->rmap->n,n=B->cmap->n; 1708 Mat Cmat; 1709 1710 PetscFunctionBegin; 1711 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); 1712 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1713 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1714 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1715 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1716 1717 *C = Cmat; 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1723 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1724 { 1725 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1726 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1727 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1728 PetscBLASInt m,n,k; 1729 PetscScalar _DOne=1.0,_DZero=0.0; 1730 PetscErrorCode ierr; 1731 1732 PetscFunctionBegin; 1733 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1734 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1735 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1736 PetscStackCallBLAS("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1737 PetscFunctionReturn(0); 1738 } 1739 1740 #undef __FUNCT__ 1741 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1742 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1743 { 1744 PetscErrorCode ierr; 1745 1746 PetscFunctionBegin; 1747 if (scall == MAT_INITIAL_MATRIX) { 1748 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1749 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1750 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 1751 } 1752 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1753 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1754 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 1755 PetscFunctionReturn(0); 1756 } 1757 1758 #undef __FUNCT__ 1759 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1760 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1761 { 1762 PetscErrorCode ierr; 1763 PetscInt m=A->cmap->n,n=B->cmap->n; 1764 Mat Cmat; 1765 1766 PetscFunctionBegin; 1767 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); 1768 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1769 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1770 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1771 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1772 1773 Cmat->assembled = PETSC_TRUE; 1774 1775 *C = Cmat; 1776 PetscFunctionReturn(0); 1777 } 1778 1779 #undef __FUNCT__ 1780 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1781 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1782 { 1783 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1784 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1785 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1786 PetscBLASInt m,n,k; 1787 PetscScalar _DOne=1.0,_DZero=0.0; 1788 PetscErrorCode ierr; 1789 1790 PetscFunctionBegin; 1791 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1792 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1793 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1794 /* 1795 Note the m and n arguments below are the number rows and columns of A', not A! 1796 */ 1797 PetscStackCallBLAS("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1798 PetscFunctionReturn(0); 1799 } 1800 1801 #undef __FUNCT__ 1802 #define __FUNCT__ "MatGetRowMax_SeqDense" 1803 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1804 { 1805 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1806 PetscErrorCode ierr; 1807 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1808 PetscScalar *x; 1809 MatScalar *aa = a->v; 1810 1811 PetscFunctionBegin; 1812 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1813 1814 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1815 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1816 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1817 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1818 for (i=0; i<m; i++) { 1819 x[i] = aa[i]; if (idx) idx[i] = 0; 1820 for (j=1; j<n; j++) { 1821 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1822 } 1823 } 1824 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1825 PetscFunctionReturn(0); 1826 } 1827 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1830 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1831 { 1832 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1833 PetscErrorCode ierr; 1834 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1835 PetscScalar *x; 1836 PetscReal atmp; 1837 MatScalar *aa = a->v; 1838 1839 PetscFunctionBegin; 1840 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1841 1842 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1843 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1844 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1845 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1846 for (i=0; i<m; i++) { 1847 x[i] = PetscAbsScalar(aa[i]); 1848 for (j=1; j<n; j++) { 1849 atmp = PetscAbsScalar(aa[i+m*j]); 1850 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1851 } 1852 } 1853 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1854 PetscFunctionReturn(0); 1855 } 1856 1857 #undef __FUNCT__ 1858 #define __FUNCT__ "MatGetRowMin_SeqDense" 1859 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1860 { 1861 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1862 PetscErrorCode ierr; 1863 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1864 PetscScalar *x; 1865 MatScalar *aa = a->v; 1866 1867 PetscFunctionBegin; 1868 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1869 1870 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1871 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1872 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1873 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1874 for (i=0; i<m; i++) { 1875 x[i] = aa[i]; if (idx) idx[i] = 0; 1876 for (j=1; j<n; j++) { 1877 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1878 } 1879 } 1880 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1881 PetscFunctionReturn(0); 1882 } 1883 1884 #undef __FUNCT__ 1885 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1886 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1887 { 1888 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1889 PetscErrorCode ierr; 1890 PetscScalar *x; 1891 1892 PetscFunctionBegin; 1893 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1894 1895 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1896 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1897 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1898 PetscFunctionReturn(0); 1899 } 1900 1901 1902 #undef __FUNCT__ 1903 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1904 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1905 { 1906 PetscErrorCode ierr; 1907 PetscInt i,j,m,n; 1908 PetscScalar *a; 1909 1910 PetscFunctionBegin; 1911 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1912 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1913 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1914 if (type == NORM_2) { 1915 for (i=0; i<n; i++) { 1916 for (j=0; j<m; j++) { 1917 norms[i] += PetscAbsScalar(a[j]*a[j]); 1918 } 1919 a += m; 1920 } 1921 } else if (type == NORM_1) { 1922 for (i=0; i<n; i++) { 1923 for (j=0; j<m; j++) { 1924 norms[i] += PetscAbsScalar(a[j]); 1925 } 1926 a += m; 1927 } 1928 } else if (type == NORM_INFINITY) { 1929 for (i=0; i<n; i++) { 1930 for (j=0; j<m; j++) { 1931 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1932 } 1933 a += m; 1934 } 1935 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1936 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1937 if (type == NORM_2) { 1938 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1939 } 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "MatSetRandom_SeqDense" 1945 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1946 { 1947 PetscErrorCode ierr; 1948 PetscScalar *a; 1949 PetscInt m,n,i; 1950 1951 PetscFunctionBegin; 1952 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1953 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1954 for (i=0; i<m*n; i++) { 1955 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1956 } 1957 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1958 PetscFunctionReturn(0); 1959 } 1960 1961 1962 /* -------------------------------------------------------------------*/ 1963 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 1964 MatGetRow_SeqDense, 1965 MatRestoreRow_SeqDense, 1966 MatMult_SeqDense, 1967 /* 4*/ MatMultAdd_SeqDense, 1968 MatMultTranspose_SeqDense, 1969 MatMultTransposeAdd_SeqDense, 1970 0, 1971 0, 1972 0, 1973 /* 10*/ 0, 1974 MatLUFactor_SeqDense, 1975 MatCholeskyFactor_SeqDense, 1976 MatSOR_SeqDense, 1977 MatTranspose_SeqDense, 1978 /* 15*/ MatGetInfo_SeqDense, 1979 MatEqual_SeqDense, 1980 MatGetDiagonal_SeqDense, 1981 MatDiagonalScale_SeqDense, 1982 MatNorm_SeqDense, 1983 /* 20*/ MatAssemblyBegin_SeqDense, 1984 MatAssemblyEnd_SeqDense, 1985 MatSetOption_SeqDense, 1986 MatZeroEntries_SeqDense, 1987 /* 24*/ MatZeroRows_SeqDense, 1988 0, 1989 0, 1990 0, 1991 0, 1992 /* 29*/ MatSetUp_SeqDense, 1993 0, 1994 0, 1995 0, 1996 0, 1997 /* 34*/ MatDuplicate_SeqDense, 1998 0, 1999 0, 2000 0, 2001 0, 2002 /* 39*/ MatAXPY_SeqDense, 2003 MatGetSubMatrices_SeqDense, 2004 0, 2005 MatGetValues_SeqDense, 2006 MatCopy_SeqDense, 2007 /* 44*/ MatGetRowMax_SeqDense, 2008 MatScale_SeqDense, 2009 0, 2010 0, 2011 0, 2012 /* 49*/ MatSetRandom_SeqDense, 2013 0, 2014 0, 2015 0, 2016 0, 2017 /* 54*/ 0, 2018 0, 2019 0, 2020 0, 2021 0, 2022 /* 59*/ 0, 2023 MatDestroy_SeqDense, 2024 MatView_SeqDense, 2025 0, 2026 0, 2027 /* 64*/ 0, 2028 0, 2029 0, 2030 0, 2031 0, 2032 /* 69*/ MatGetRowMaxAbs_SeqDense, 2033 0, 2034 0, 2035 0, 2036 0, 2037 /* 74*/ 0, 2038 0, 2039 0, 2040 0, 2041 0, 2042 /* 79*/ 0, 2043 0, 2044 0, 2045 0, 2046 /* 83*/ MatLoad_SeqDense, 2047 0, 2048 MatIsHermitian_SeqDense, 2049 0, 2050 0, 2051 0, 2052 /* 89*/ MatMatMult_SeqDense_SeqDense, 2053 MatMatMultSymbolic_SeqDense_SeqDense, 2054 MatMatMultNumeric_SeqDense_SeqDense, 2055 0, 2056 0, 2057 /* 94*/ 0, 2058 0, 2059 0, 2060 0, 2061 0, 2062 /* 99*/ 0, 2063 0, 2064 0, 2065 MatConjugate_SeqDense, 2066 0, 2067 /*104*/ 0, 2068 MatRealPart_SeqDense, 2069 MatImaginaryPart_SeqDense, 2070 0, 2071 0, 2072 /*109*/ MatMatSolve_SeqDense, 2073 0, 2074 MatGetRowMin_SeqDense, 2075 MatGetColumnVector_SeqDense, 2076 0, 2077 /*114*/ 0, 2078 0, 2079 0, 2080 0, 2081 0, 2082 /*119*/ 0, 2083 0, 2084 0, 2085 0, 2086 0, 2087 /*124*/ 0, 2088 MatGetColumnNorms_SeqDense, 2089 0, 2090 0, 2091 0, 2092 /*129*/ 0, 2093 MatTransposeMatMult_SeqDense_SeqDense, 2094 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2095 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2096 0, 2097 /*134*/ 0, 2098 0, 2099 0, 2100 0, 2101 0, 2102 /*139*/ 0, 2103 0, 2104 0 2105 }; 2106 2107 #undef __FUNCT__ 2108 #define __FUNCT__ "MatCreateSeqDense" 2109 /*@C 2110 MatCreateSeqDense - Creates a sequential dense matrix that 2111 is stored in column major order (the usual Fortran 77 manner). Many 2112 of the matrix operations use the BLAS and LAPACK routines. 2113 2114 Collective on MPI_Comm 2115 2116 Input Parameters: 2117 + comm - MPI communicator, set to PETSC_COMM_SELF 2118 . m - number of rows 2119 . n - number of columns 2120 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2121 to control all matrix memory allocation. 2122 2123 Output Parameter: 2124 . A - the matrix 2125 2126 Notes: 2127 The data input variable is intended primarily for Fortran programmers 2128 who wish to allocate their own matrix memory space. Most users should 2129 set data=NULL. 2130 2131 Level: intermediate 2132 2133 .keywords: dense, matrix, LAPACK, BLAS 2134 2135 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2136 @*/ 2137 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2138 { 2139 PetscErrorCode ierr; 2140 2141 PetscFunctionBegin; 2142 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2143 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2144 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2145 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2146 PetscFunctionReturn(0); 2147 } 2148 2149 #undef __FUNCT__ 2150 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2151 /*@C 2152 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2153 2154 Collective on MPI_Comm 2155 2156 Input Parameters: 2157 + A - the matrix 2158 - data - the array (or NULL) 2159 2160 Notes: 2161 The data input variable is intended primarily for Fortran programmers 2162 who wish to allocate their own matrix memory space. Most users should 2163 need not call this routine. 2164 2165 Level: intermediate 2166 2167 .keywords: dense, matrix, LAPACK, BLAS 2168 2169 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2170 2171 @*/ 2172 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2173 { 2174 PetscErrorCode ierr; 2175 2176 PetscFunctionBegin; 2177 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2178 PetscFunctionReturn(0); 2179 } 2180 2181 #undef __FUNCT__ 2182 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2183 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2184 { 2185 Mat_SeqDense *b; 2186 PetscErrorCode ierr; 2187 2188 PetscFunctionBegin; 2189 B->preallocated = PETSC_TRUE; 2190 2191 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2192 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2193 2194 b = (Mat_SeqDense*)B->data; 2195 b->Mmax = B->rmap->n; 2196 b->Nmax = B->cmap->n; 2197 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2198 2199 if (!data) { /* petsc-allocated storage */ 2200 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2201 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2202 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2203 2204 b->user_alloc = PETSC_FALSE; 2205 } else { /* user-allocated storage */ 2206 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2207 b->v = data; 2208 b->user_alloc = PETSC_TRUE; 2209 } 2210 B->assembled = PETSC_TRUE; 2211 PetscFunctionReturn(0); 2212 } 2213 2214 #undef __FUNCT__ 2215 #define __FUNCT__ "MatSeqDenseSetLDA" 2216 /*@C 2217 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2218 2219 Input parameter: 2220 + A - the matrix 2221 - lda - the leading dimension 2222 2223 Notes: 2224 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2225 it asserts that the preallocation has a leading dimension (the LDA parameter 2226 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2227 2228 Level: intermediate 2229 2230 .keywords: dense, matrix, LAPACK, BLAS 2231 2232 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2233 2234 @*/ 2235 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2236 { 2237 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2238 2239 PetscFunctionBegin; 2240 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); 2241 b->lda = lda; 2242 b->changelda = PETSC_FALSE; 2243 b->Mmax = PetscMax(b->Mmax,lda); 2244 PetscFunctionReturn(0); 2245 } 2246 2247 /*MC 2248 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2249 2250 Options Database Keys: 2251 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2252 2253 Level: beginner 2254 2255 .seealso: MatCreateSeqDense() 2256 2257 M*/ 2258 2259 #undef __FUNCT__ 2260 #define __FUNCT__ "MatCreate_SeqDense" 2261 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2262 { 2263 Mat_SeqDense *b; 2264 PetscErrorCode ierr; 2265 PetscMPIInt size; 2266 2267 PetscFunctionBegin; 2268 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2269 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2270 2271 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2272 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2273 B->data = (void*)b; 2274 2275 b->pivots = 0; 2276 b->roworiented = PETSC_TRUE; 2277 b->v = 0; 2278 b->changelda = PETSC_FALSE; 2279 2280 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2281 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2282 2283 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2284 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2285 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2288 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2289 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2292 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2293 PetscFunctionReturn(0); 2294 } 2295