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 EXTERN_C_BEGIN 11 #undef __FUNCT__ 12 #define __FUNCT__ "MatConvert_SeqDense_SeqAIJ" 13 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; 19 PetscInt *rows; 20 MatScalar *aa = a->v; 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 = MatSeqAIJSetPreallocation(B,A->cmap->n,NULL);CHKERRQ(ierr); 27 28 ierr = PetscMalloc(A->rmap->n*sizeof(PetscInt),&rows);CHKERRQ(ierr); 29 for (i=0; i<A->rmap->n; i++) rows[i] = i; 30 31 for (i=0; i<A->cmap->n; i++) { 32 ierr = MatSetValues(B,A->rmap->n,rows,1,&i,aa,INSERT_VALUES);CHKERRQ(ierr); 33 aa += a->lda; 34 } 35 ierr = PetscFree(rows);CHKERRQ(ierr); 36 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 37 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 38 39 if (reuse == MAT_REUSE_MATRIX) { 40 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 41 } else { 42 *newmat = B; 43 } 44 PetscFunctionReturn(0); 45 } 46 EXTERN_C_END 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatAXPY_SeqDense" 50 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 51 { 52 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 53 PetscScalar oalpha = alpha; 54 PetscInt j; 55 PetscBLASInt N,m,ldax,lday,one = 1; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 60 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 61 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 62 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 63 if (ldax>m || lday>m) { 64 for (j=0; j<X->cmap->n; j++) { 65 PetscStackCall("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one)); 66 } 67 } else { 68 PetscStackCall("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one)); 69 } 70 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "MatGetInfo_SeqDense" 76 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 77 { 78 PetscInt N = A->rmap->n*A->cmap->n; 79 80 PetscFunctionBegin; 81 info->block_size = 1.0; 82 info->nz_allocated = (double)N; 83 info->nz_used = (double)N; 84 info->nz_unneeded = (double)0; 85 info->assemblies = (double)A->num_ass; 86 info->mallocs = 0; 87 info->memory = ((PetscObject)A)->mem; 88 info->fill_ratio_given = 0; 89 info->fill_ratio_needed = 0; 90 info->factor_mallocs = 0; 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatScale_SeqDense" 96 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 97 { 98 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 99 PetscScalar oalpha = alpha; 100 PetscErrorCode ierr; 101 PetscBLASInt one = 1,j,nz,lda; 102 103 PetscFunctionBegin; 104 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 105 if (lda>A->rmap->n) { 106 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 107 for (j=0; j<A->cmap->n; j++) { 108 PetscStackCall("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one)); 109 } 110 } else { 111 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 112 PetscStackCall("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one)); 113 } 114 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "MatIsHermitian_SeqDense" 120 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 121 { 122 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 123 PetscInt i,j,m = A->rmap->n,N; 124 PetscScalar *v = a->v; 125 126 PetscFunctionBegin; 127 *fl = PETSC_FALSE; 128 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 129 N = a->lda; 130 131 for (i=0; i<m; i++) { 132 for (j=i+1; j<m; j++) { 133 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 134 } 135 } 136 *fl = PETSC_TRUE; 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 142 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 143 { 144 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 145 PetscErrorCode ierr; 146 PetscInt lda = (PetscInt)mat->lda,j,m; 147 148 PetscFunctionBegin; 149 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 150 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 151 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 152 if (cpvalues == MAT_COPY_VALUES) { 153 l = (Mat_SeqDense*)newi->data; 154 if (lda>A->rmap->n) { 155 m = A->rmap->n; 156 for (j=0; j<A->cmap->n; j++) { 157 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 158 } 159 } else { 160 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 161 } 162 } 163 newi->assembled = PETSC_TRUE; 164 PetscFunctionReturn(0); 165 } 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "MatDuplicate_SeqDense" 169 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 170 { 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 175 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 176 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 177 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 182 extern PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 186 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 187 { 188 MatFactorInfo info; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 193 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "MatSolve_SeqDense" 199 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 200 { 201 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 202 PetscErrorCode ierr; 203 PetscScalar *x,*y; 204 PetscBLASInt one = 1,info,m; 205 206 PetscFunctionBegin; 207 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 208 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 209 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 210 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 211 if (A->factortype == MAT_FACTOR_LU) { 212 #if defined(PETSC_MISSING_LAPACK_GETRS) 213 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 214 #else 215 PetscStackCall("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 216 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 217 #endif 218 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 219 #if defined(PETSC_MISSING_LAPACK_POTRS) 220 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 221 #else 222 PetscStackCall("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 223 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 224 #endif 225 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 226 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 227 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 228 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatMatSolve_SeqDense" 234 PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 235 { 236 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 237 PetscErrorCode ierr; 238 PetscScalar *b,*x; 239 PetscInt n; 240 PetscBLASInt nrhs,info,m; 241 PetscBool flg; 242 243 PetscFunctionBegin; 244 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 245 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 246 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 247 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 248 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 249 250 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 251 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 252 ierr = MatDenseGetArray(B,&b);CHKERRQ(ierr); 253 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 254 255 ierr = PetscMemcpy(x,b,m*nrhs*sizeof(PetscScalar));CHKERRQ(ierr); 256 257 if (A->factortype == MAT_FACTOR_LU) { 258 #if defined(PETSC_MISSING_LAPACK_GETRS) 259 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 260 #else 261 PetscStackCall("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 262 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 263 #endif 264 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 265 #if defined(PETSC_MISSING_LAPACK_POTRS) 266 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 267 #else 268 PetscStackCall("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 269 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 270 #endif 271 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 272 273 ierr = MatDenseRestoreArray(B,&b);CHKERRQ(ierr); 274 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 275 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "MatSolveTranspose_SeqDense" 281 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 282 { 283 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 284 PetscErrorCode ierr; 285 PetscScalar *x,*y; 286 PetscBLASInt one = 1,info,m; 287 288 PetscFunctionBegin; 289 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 290 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 291 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 292 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 293 /* assume if pivots exist then use LU; else Cholesky */ 294 if (mat->pivots) { 295 #if defined(PETSC_MISSING_LAPACK_GETRS) 296 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 297 #else 298 PetscStackCall("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 299 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 300 #endif 301 } else { 302 #if defined(PETSC_MISSING_LAPACK_POTRS) 303 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 304 #else 305 PetscStackCall("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 306 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 307 #endif 308 } 309 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 310 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 311 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "MatSolveAdd_SeqDense" 317 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 318 { 319 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 320 PetscErrorCode ierr; 321 PetscScalar *x,*y,sone = 1.0; 322 Vec tmp = 0; 323 PetscBLASInt one = 1,info,m; 324 325 PetscFunctionBegin; 326 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 327 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 328 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 329 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 330 if (yy == zz) { 331 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 332 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 333 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 334 } 335 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 336 /* assume if pivots exist then use LU; else Cholesky */ 337 if (mat->pivots) { 338 #if defined(PETSC_MISSING_LAPACK_GETRS) 339 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 340 #else 341 PetscStackCall("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 342 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 343 #endif 344 } else { 345 #if defined(PETSC_MISSING_LAPACK_POTRS) 346 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 347 #else 348 PetscStackCall("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 349 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 350 #endif 351 } 352 if (tmp) { 353 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 354 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 355 } else { 356 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 357 } 358 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 359 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 360 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 366 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 367 { 368 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 369 PetscErrorCode ierr; 370 PetscScalar *x,*y,sone = 1.0; 371 Vec tmp; 372 PetscBLASInt one = 1,info,m; 373 374 PetscFunctionBegin; 375 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 376 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 377 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 378 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 379 if (yy == zz) { 380 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 381 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 382 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 383 } 384 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 385 /* assume if pivots exist then use LU; else Cholesky */ 386 if (mat->pivots) { 387 #if defined(PETSC_MISSING_LAPACK_GETRS) 388 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 389 #else 390 PetscStackCall("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 391 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 392 #endif 393 } else { 394 #if defined(PETSC_MISSING_LAPACK_POTRS) 395 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 396 #else 397 PetscStackCall("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 398 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad solve"); 399 #endif 400 } 401 if (tmp) { 402 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 403 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 404 } else { 405 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 406 } 407 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 408 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 409 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 /* ---------------------------------------------------------------*/ 414 /* COMMENT: I have chosen to hide row permutation in the pivots, 415 rather than put it in the Mat->row slot.*/ 416 #undef __FUNCT__ 417 #define __FUNCT__ "MatLUFactor_SeqDense" 418 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 419 { 420 #if defined(PETSC_MISSING_LAPACK_GETRF) 421 PetscFunctionBegin; 422 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 423 #else 424 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 425 PetscErrorCode ierr; 426 PetscBLASInt n,m,info; 427 428 PetscFunctionBegin; 429 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 430 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 431 if (!mat->pivots) { 432 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 433 ierr = PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 434 } 435 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 436 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 437 PetscStackCall("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 438 ierr = PetscFPTrapPop();CHKERRQ(ierr); 439 440 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 441 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 442 A->ops->solve = MatSolve_SeqDense; 443 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 444 A->ops->solveadd = MatSolveAdd_SeqDense; 445 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 446 A->factortype = MAT_FACTOR_LU; 447 448 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 449 #endif 450 PetscFunctionReturn(0); 451 } 452 453 #undef __FUNCT__ 454 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 455 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 456 { 457 #if defined(PETSC_MISSING_LAPACK_POTRF) 458 PetscFunctionBegin; 459 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 460 #else 461 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 462 PetscErrorCode ierr; 463 PetscBLASInt info,n; 464 465 PetscFunctionBegin; 466 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 467 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 468 469 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 470 PetscStackCall("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 471 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 472 A->ops->solve = MatSolve_SeqDense; 473 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 474 A->ops->solveadd = MatSolveAdd_SeqDense; 475 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 476 A->factortype = MAT_FACTOR_CHOLESKY; 477 478 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 479 #endif 480 PetscFunctionReturn(0); 481 } 482 483 484 #undef __FUNCT__ 485 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 486 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 487 { 488 PetscErrorCode ierr; 489 MatFactorInfo info; 490 491 PetscFunctionBegin; 492 info.fill = 1.0; 493 494 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 495 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 #undef __FUNCT__ 500 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 501 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 502 { 503 PetscFunctionBegin; 504 fact->assembled = PETSC_TRUE; 505 fact->preallocated = PETSC_TRUE; 506 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 507 PetscFunctionReturn(0); 508 } 509 510 #undef __FUNCT__ 511 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 512 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 513 { 514 PetscFunctionBegin; 515 fact->preallocated = PETSC_TRUE; 516 fact->assembled = PETSC_TRUE; 517 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 518 PetscFunctionReturn(0); 519 } 520 521 EXTERN_C_BEGIN 522 #undef __FUNCT__ 523 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 524 PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 525 { 526 PetscErrorCode ierr; 527 528 PetscFunctionBegin; 529 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 530 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 531 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 532 if (ftype == MAT_FACTOR_LU) { 533 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 534 } else { 535 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 536 } 537 (*fact)->factortype = ftype; 538 PetscFunctionReturn(0); 539 } 540 EXTERN_C_END 541 542 /* ------------------------------------------------------------------*/ 543 #undef __FUNCT__ 544 #define __FUNCT__ "MatSOR_SeqDense" 545 PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 546 { 547 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 548 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 549 PetscErrorCode ierr; 550 PetscInt m = A->rmap->n,i; 551 PetscBLASInt o = 1,bm; 552 553 PetscFunctionBegin; 554 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 555 if (flag & SOR_ZERO_INITIAL_GUESS) { 556 /* this is a hack fix, should have another version without the second BLASdot */ 557 ierr = VecSet(xx,zero);CHKERRQ(ierr); 558 } 559 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 560 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 561 its = its*lits; 562 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 563 while (its--) { 564 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 565 for (i=0; i<m; i++) { 566 PetscStackCall("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 567 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 568 } 569 } 570 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 571 for (i=m-1; i>=0; i--) { 572 PetscStackCall("BLASdot",xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o)); 573 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 574 } 575 } 576 } 577 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 578 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 579 PetscFunctionReturn(0); 580 } 581 582 /* -----------------------------------------------------------------*/ 583 #undef __FUNCT__ 584 #define __FUNCT__ "MatMultTranspose_SeqDense" 585 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 586 { 587 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 588 PetscScalar *v = mat->v,*x,*y; 589 PetscErrorCode ierr; 590 PetscBLASInt m, n,_One=1; 591 PetscScalar _DOne=1.0,_DZero=0.0; 592 593 PetscFunctionBegin; 594 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 595 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 596 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 597 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 598 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 599 PetscStackCall("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 600 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 601 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 602 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "MatMult_SeqDense" 608 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 609 { 610 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 611 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 612 PetscErrorCode ierr; 613 PetscBLASInt m, n, _One=1; 614 615 PetscFunctionBegin; 616 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 617 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 618 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 619 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 620 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 621 PetscStackCall("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 622 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 623 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 624 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 625 PetscFunctionReturn(0); 626 } 627 628 #undef __FUNCT__ 629 #define __FUNCT__ "MatMultAdd_SeqDense" 630 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 631 { 632 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 633 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 634 PetscErrorCode ierr; 635 PetscBLASInt m, n, _One=1; 636 637 PetscFunctionBegin; 638 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 639 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 640 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 641 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 642 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 643 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 644 PetscStackCall("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 645 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 646 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 647 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 648 PetscFunctionReturn(0); 649 } 650 651 #undef __FUNCT__ 652 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 653 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 654 { 655 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 656 PetscScalar *v = mat->v,*x,*y; 657 PetscErrorCode ierr; 658 PetscBLASInt m, n, _One=1; 659 PetscScalar _DOne=1.0; 660 661 PetscFunctionBegin; 662 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 663 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 664 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 665 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 666 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 667 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 668 PetscStackCall("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 669 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 670 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 671 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 675 /* -----------------------------------------------------------------*/ 676 #undef __FUNCT__ 677 #define __FUNCT__ "MatGetRow_SeqDense" 678 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 679 { 680 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 681 PetscScalar *v; 682 PetscErrorCode ierr; 683 PetscInt i; 684 685 PetscFunctionBegin; 686 *ncols = A->cmap->n; 687 if (cols) { 688 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);CHKERRQ(ierr); 689 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 690 } 691 if (vals) { 692 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 693 v = mat->v + row; 694 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 695 } 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "MatRestoreRow_SeqDense" 701 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 702 { 703 PetscErrorCode ierr; 704 705 PetscFunctionBegin; 706 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 707 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 708 PetscFunctionReturn(0); 709 } 710 /* ----------------------------------------------------------------*/ 711 #undef __FUNCT__ 712 #define __FUNCT__ "MatSetValues_SeqDense" 713 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 714 { 715 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 716 PetscInt i,j,idx=0; 717 718 PetscFunctionBegin; 719 if (v) PetscValidScalarPointer(v,6); 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 = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&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 = PetscMalloc((M+1)*sizeof(PetscInt),&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 = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 864 cols = scols; 865 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 866 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&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 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");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,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 912 } else if (PetscRealPart(*v)) { 913 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 914 } 915 #else 916 if (*v) { 917 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*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 } else { 940 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 941 } 942 943 for (i=0; i<A->rmap->n; i++) { 944 v = a->v + i; 945 for (j=0; j<A->cmap->n; j++) { 946 #if defined(PETSC_USE_COMPLEX) 947 if (allreal) { 948 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 949 } else { 950 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 951 } 952 #else 953 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 954 #endif 955 v += a->lda; 956 } 957 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 958 } 959 if (format == PETSC_VIEWER_ASCII_MATLAB) { 960 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 961 } 962 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 963 } 964 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 965 PetscFunctionReturn(0); 966 } 967 968 #undef __FUNCT__ 969 #define __FUNCT__ "MatView_SeqDense_Binary" 970 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 971 { 972 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 973 PetscErrorCode ierr; 974 int fd; 975 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 976 PetscScalar *v,*anonz,*vals; 977 PetscViewerFormat format; 978 979 PetscFunctionBegin; 980 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 981 982 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 983 if (format == PETSC_VIEWER_NATIVE) { 984 /* store the matrix as a dense matrix */ 985 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 986 987 col_lens[0] = MAT_FILE_CLASSID; 988 col_lens[1] = m; 989 col_lens[2] = n; 990 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 991 992 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 993 ierr = PetscFree(col_lens);CHKERRQ(ierr); 994 995 /* write out matrix, by rows */ 996 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 997 v = a->v; 998 for (j=0; j<n; j++) { 999 for (i=0; i<m; i++) { 1000 vals[j + i*n] = *v++; 1001 } 1002 } 1003 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1004 ierr = PetscFree(vals);CHKERRQ(ierr); 1005 } else { 1006 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 1007 1008 col_lens[0] = MAT_FILE_CLASSID; 1009 col_lens[1] = m; 1010 col_lens[2] = n; 1011 col_lens[3] = nz; 1012 1013 /* store lengths of each row and write (including header) to file */ 1014 for (i=0; i<m; i++) col_lens[4+i] = n; 1015 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1016 1017 /* Possibly should write in smaller increments, not whole matrix at once? */ 1018 /* store column indices (zero start index) */ 1019 ict = 0; 1020 for (i=0; i<m; i++) { 1021 for (j=0; j<n; j++) col_lens[ict++] = j; 1022 } 1023 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1024 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1025 1026 /* store nonzero values */ 1027 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 1028 ict = 0; 1029 for (i=0; i<m; i++) { 1030 v = a->v + i; 1031 for (j=0; j<n; j++) { 1032 anonz[ict++] = *v; v += a->lda; 1033 } 1034 } 1035 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1036 ierr = PetscFree(anonz);CHKERRQ(ierr); 1037 } 1038 PetscFunctionReturn(0); 1039 } 1040 1041 #include <petscdraw.h> 1042 #undef __FUNCT__ 1043 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 1044 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1045 { 1046 Mat A = (Mat) Aa; 1047 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1048 PetscErrorCode ierr; 1049 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 1050 PetscScalar *v = a->v; 1051 PetscViewer viewer; 1052 PetscDraw popup; 1053 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 1054 PetscViewerFormat format; 1055 1056 PetscFunctionBegin; 1057 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1058 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1059 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1060 1061 /* Loop over matrix elements drawing boxes */ 1062 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1063 /* Blue for negative and Red for positive */ 1064 color = PETSC_DRAW_BLUE; 1065 for (j = 0; j < n; j++) { 1066 x_l = j; 1067 x_r = x_l + 1.0; 1068 for (i = 0; i < m; i++) { 1069 y_l = m - i - 1.0; 1070 y_r = y_l + 1.0; 1071 if (PetscRealPart(v[j*m+i]) > 0.) { 1072 color = PETSC_DRAW_RED; 1073 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1074 color = PETSC_DRAW_BLUE; 1075 } else { 1076 continue; 1077 } 1078 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1079 } 1080 } 1081 } else { 1082 /* use contour shading to indicate magnitude of values */ 1083 /* first determine max of all nonzero values */ 1084 for (i = 0; i < m*n; i++) { 1085 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1086 } 1087 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1088 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1089 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1090 for (j = 0; j < n; j++) { 1091 x_l = j; 1092 x_r = x_l + 1.0; 1093 for (i = 0; i < m; i++) { 1094 y_l = m - i - 1.0; 1095 y_r = y_l + 1.0; 1096 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1097 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1098 } 1099 } 1100 } 1101 PetscFunctionReturn(0); 1102 } 1103 1104 #undef __FUNCT__ 1105 #define __FUNCT__ "MatView_SeqDense_Draw" 1106 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1107 { 1108 PetscDraw draw; 1109 PetscBool isnull; 1110 PetscReal xr,yr,xl,yl,h,w; 1111 PetscErrorCode ierr; 1112 1113 PetscFunctionBegin; 1114 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1115 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1116 if (isnull) PetscFunctionReturn(0); 1117 1118 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1119 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1120 xr += w; yr += h; xl = -w; yl = -h; 1121 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1122 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1123 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 #undef __FUNCT__ 1128 #define __FUNCT__ "MatView_SeqDense" 1129 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1130 { 1131 PetscErrorCode ierr; 1132 PetscBool iascii,isbinary,isdraw; 1133 1134 PetscFunctionBegin; 1135 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1136 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1137 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1138 1139 if (iascii) { 1140 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1141 } else if (isbinary) { 1142 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1143 } else if (isdraw) { 1144 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "MatDestroy_SeqDense" 1151 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1152 { 1153 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1154 PetscErrorCode ierr; 1155 1156 PetscFunctionBegin; 1157 #if defined(PETSC_USE_LOG) 1158 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1159 #endif 1160 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1161 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1162 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1163 1164 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1165 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseGetArray_C","",NULL);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatDenseRestoreArray_C","",NULL);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",NULL);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",NULL);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",NULL);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_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 = PetscMalloc((n+1)*sizeof(Mat),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 = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1693 } 1694 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1700 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1701 { 1702 PetscErrorCode ierr; 1703 PetscInt m=A->rmap->n,n=B->cmap->n; 1704 Mat Cmat; 1705 1706 PetscFunctionBegin; 1707 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); 1708 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1709 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1710 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1711 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1712 1713 *C = Cmat; 1714 PetscFunctionReturn(0); 1715 } 1716 1717 #undef __FUNCT__ 1718 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1719 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1720 { 1721 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1722 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1723 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1724 PetscBLASInt m,n,k; 1725 PetscScalar _DOne=1.0,_DZero=0.0; 1726 PetscErrorCode ierr; 1727 1728 PetscFunctionBegin; 1729 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 1730 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1731 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 1732 PetscStackCall("BLASgemv",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1733 PetscFunctionReturn(0); 1734 } 1735 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "MatTransposeMatMult_SeqDense_SeqDense" 1738 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1739 { 1740 PetscErrorCode ierr; 1741 1742 PetscFunctionBegin; 1743 if (scall == MAT_INITIAL_MATRIX) { 1744 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1745 } 1746 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1747 PetscFunctionReturn(0); 1748 } 1749 1750 #undef __FUNCT__ 1751 #define __FUNCT__ "MatTransposeMatMultSymbolic_SeqDense_SeqDense" 1752 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1753 { 1754 PetscErrorCode ierr; 1755 PetscInt m=A->cmap->n,n=B->cmap->n; 1756 Mat Cmat; 1757 1758 PetscFunctionBegin; 1759 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); 1760 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1761 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1762 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1763 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 1764 1765 Cmat->assembled = PETSC_TRUE; 1766 1767 *C = Cmat; 1768 PetscFunctionReturn(0); 1769 } 1770 1771 #undef __FUNCT__ 1772 #define __FUNCT__ "MatTransposeMatMultNumeric_SeqDense_SeqDense" 1773 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1774 { 1775 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1776 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1777 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1778 PetscBLASInt m,n,k; 1779 PetscScalar _DOne=1.0,_DZero=0.0; 1780 PetscErrorCode ierr; 1781 1782 PetscFunctionBegin; 1783 ierr = PetscBLASIntCast(A->cmap->n,&m);CHKERRQ(ierr); 1784 ierr = PetscBLASIntCast(B->cmap->n,&n);CHKERRQ(ierr); 1785 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 1786 /* 1787 Note the m and n arguments below are the number rows and columns of A', not A! 1788 */ 1789 PetscStackCall("BLASgemv",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 1790 PetscFunctionReturn(0); 1791 } 1792 1793 #undef __FUNCT__ 1794 #define __FUNCT__ "MatGetRowMax_SeqDense" 1795 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1796 { 1797 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1798 PetscErrorCode ierr; 1799 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1800 PetscScalar *x; 1801 MatScalar *aa = a->v; 1802 1803 PetscFunctionBegin; 1804 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1805 1806 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1807 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1808 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1809 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1810 for (i=0; i<m; i++) { 1811 x[i] = aa[i]; if (idx) idx[i] = 0; 1812 for (j=1; j<n; j++) { 1813 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1814 } 1815 } 1816 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 #undef __FUNCT__ 1821 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1822 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1823 { 1824 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1825 PetscErrorCode ierr; 1826 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1827 PetscScalar *x; 1828 PetscReal atmp; 1829 MatScalar *aa = a->v; 1830 1831 PetscFunctionBegin; 1832 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1833 1834 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1835 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1836 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1837 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1838 for (i=0; i<m; i++) { 1839 x[i] = PetscAbsScalar(aa[i]); 1840 for (j=1; j<n; j++) { 1841 atmp = PetscAbsScalar(aa[i+m*j]); 1842 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1843 } 1844 } 1845 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1846 PetscFunctionReturn(0); 1847 } 1848 1849 #undef __FUNCT__ 1850 #define __FUNCT__ "MatGetRowMin_SeqDense" 1851 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1852 { 1853 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1854 PetscErrorCode ierr; 1855 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1856 PetscScalar *x; 1857 MatScalar *aa = a->v; 1858 1859 PetscFunctionBegin; 1860 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1861 1862 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1863 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1864 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1865 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1866 for (i=0; i<m; i++) { 1867 x[i] = aa[i]; if (idx) idx[i] = 0; 1868 for (j=1; j<n; j++) { 1869 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1870 } 1871 } 1872 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1878 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1879 { 1880 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1881 PetscErrorCode ierr; 1882 PetscScalar *x; 1883 1884 PetscFunctionBegin; 1885 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1886 1887 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1888 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1889 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 1893 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 1896 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 1897 { 1898 PetscErrorCode ierr; 1899 PetscInt i,j,m,n; 1900 PetscScalar *a; 1901 1902 PetscFunctionBegin; 1903 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 1904 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 1905 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 1906 if (type == NORM_2) { 1907 for (i=0; i<n; i++) { 1908 for (j=0; j<m; j++) { 1909 norms[i] += PetscAbsScalar(a[j]*a[j]); 1910 } 1911 a += m; 1912 } 1913 } else if (type == NORM_1) { 1914 for (i=0; i<n; i++) { 1915 for (j=0; j<m; j++) { 1916 norms[i] += PetscAbsScalar(a[j]); 1917 } 1918 a += m; 1919 } 1920 } else if (type == NORM_INFINITY) { 1921 for (i=0; i<n; i++) { 1922 for (j=0; j<m; j++) { 1923 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 1924 } 1925 a += m; 1926 } 1927 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 1928 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 1929 if (type == NORM_2) { 1930 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1931 } 1932 PetscFunctionReturn(0); 1933 } 1934 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "MatSetRandom_SeqDense" 1937 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 1938 { 1939 PetscErrorCode ierr; 1940 PetscScalar *a; 1941 PetscInt m,n,i; 1942 1943 PetscFunctionBegin; 1944 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 1945 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 1946 for (i=0; i<m*n; i++) { 1947 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 1948 } 1949 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 1954 /* -------------------------------------------------------------------*/ 1955 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 1956 MatGetRow_SeqDense, 1957 MatRestoreRow_SeqDense, 1958 MatMult_SeqDense, 1959 /* 4*/ MatMultAdd_SeqDense, 1960 MatMultTranspose_SeqDense, 1961 MatMultTransposeAdd_SeqDense, 1962 0, 1963 0, 1964 0, 1965 /* 10*/ 0, 1966 MatLUFactor_SeqDense, 1967 MatCholeskyFactor_SeqDense, 1968 MatSOR_SeqDense, 1969 MatTranspose_SeqDense, 1970 /* 15*/ MatGetInfo_SeqDense, 1971 MatEqual_SeqDense, 1972 MatGetDiagonal_SeqDense, 1973 MatDiagonalScale_SeqDense, 1974 MatNorm_SeqDense, 1975 /* 20*/ MatAssemblyBegin_SeqDense, 1976 MatAssemblyEnd_SeqDense, 1977 MatSetOption_SeqDense, 1978 MatZeroEntries_SeqDense, 1979 /* 24*/ MatZeroRows_SeqDense, 1980 0, 1981 0, 1982 0, 1983 0, 1984 /* 29*/ MatSetUp_SeqDense, 1985 0, 1986 0, 1987 0, 1988 0, 1989 /* 34*/ MatDuplicate_SeqDense, 1990 0, 1991 0, 1992 0, 1993 0, 1994 /* 39*/ MatAXPY_SeqDense, 1995 MatGetSubMatrices_SeqDense, 1996 0, 1997 MatGetValues_SeqDense, 1998 MatCopy_SeqDense, 1999 /* 44*/ MatGetRowMax_SeqDense, 2000 MatScale_SeqDense, 2001 0, 2002 0, 2003 0, 2004 /* 49*/ MatSetRandom_SeqDense, 2005 0, 2006 0, 2007 0, 2008 0, 2009 /* 54*/ 0, 2010 0, 2011 0, 2012 0, 2013 0, 2014 /* 59*/ 0, 2015 MatDestroy_SeqDense, 2016 MatView_SeqDense, 2017 0, 2018 0, 2019 /* 64*/ 0, 2020 0, 2021 0, 2022 0, 2023 0, 2024 /* 69*/ MatGetRowMaxAbs_SeqDense, 2025 0, 2026 0, 2027 0, 2028 0, 2029 /* 74*/ 0, 2030 0, 2031 0, 2032 0, 2033 0, 2034 /* 79*/ 0, 2035 0, 2036 0, 2037 0, 2038 /* 83*/ MatLoad_SeqDense, 2039 0, 2040 MatIsHermitian_SeqDense, 2041 0, 2042 0, 2043 0, 2044 /* 89*/ MatMatMult_SeqDense_SeqDense, 2045 MatMatMultSymbolic_SeqDense_SeqDense, 2046 MatMatMultNumeric_SeqDense_SeqDense, 2047 0, 2048 0, 2049 /* 94*/ 0, 2050 0, 2051 0, 2052 0, 2053 0, 2054 /* 99*/ 0, 2055 0, 2056 0, 2057 MatConjugate_SeqDense, 2058 0, 2059 /*104*/ 0, 2060 MatRealPart_SeqDense, 2061 MatImaginaryPart_SeqDense, 2062 0, 2063 0, 2064 /*109*/ MatMatSolve_SeqDense, 2065 0, 2066 MatGetRowMin_SeqDense, 2067 MatGetColumnVector_SeqDense, 2068 0, 2069 /*114*/ 0, 2070 0, 2071 0, 2072 0, 2073 0, 2074 /*119*/ 0, 2075 0, 2076 0, 2077 0, 2078 0, 2079 /*124*/ 0, 2080 MatGetColumnNorms_SeqDense, 2081 0, 2082 0, 2083 0, 2084 /*129*/ 0, 2085 MatTransposeMatMult_SeqDense_SeqDense, 2086 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2087 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2088 }; 2089 2090 #undef __FUNCT__ 2091 #define __FUNCT__ "MatCreateSeqDense" 2092 /*@C 2093 MatCreateSeqDense - Creates a sequential dense matrix that 2094 is stored in column major order (the usual Fortran 77 manner). Many 2095 of the matrix operations use the BLAS and LAPACK routines. 2096 2097 Collective on MPI_Comm 2098 2099 Input Parameters: 2100 + comm - MPI communicator, set to PETSC_COMM_SELF 2101 . m - number of rows 2102 . n - number of columns 2103 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2104 to control all matrix memory allocation. 2105 2106 Output Parameter: 2107 . A - the matrix 2108 2109 Notes: 2110 The data input variable is intended primarily for Fortran programmers 2111 who wish to allocate their own matrix memory space. Most users should 2112 set data=NULL. 2113 2114 Level: intermediate 2115 2116 .keywords: dense, matrix, LAPACK, BLAS 2117 2118 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2119 @*/ 2120 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2121 { 2122 PetscErrorCode ierr; 2123 2124 PetscFunctionBegin; 2125 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2126 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2127 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2128 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2129 PetscFunctionReturn(0); 2130 } 2131 2132 #undef __FUNCT__ 2133 #define __FUNCT__ "MatSeqDenseSetPreallocation" 2134 /*@C 2135 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2136 2137 Collective on MPI_Comm 2138 2139 Input Parameters: 2140 + A - the matrix 2141 - data - the array (or NULL) 2142 2143 Notes: 2144 The data input variable is intended primarily for Fortran programmers 2145 who wish to allocate their own matrix memory space. Most users should 2146 need not call this routine. 2147 2148 Level: intermediate 2149 2150 .keywords: dense, matrix, LAPACK, BLAS 2151 2152 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2153 2154 @*/ 2155 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2156 { 2157 PetscErrorCode ierr; 2158 2159 PetscFunctionBegin; 2160 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2161 PetscFunctionReturn(0); 2162 } 2163 2164 EXTERN_C_BEGIN 2165 #undef __FUNCT__ 2166 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 2167 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2168 { 2169 Mat_SeqDense *b; 2170 PetscErrorCode ierr; 2171 2172 PetscFunctionBegin; 2173 B->preallocated = PETSC_TRUE; 2174 2175 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2176 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2177 2178 b = (Mat_SeqDense*)B->data; 2179 b->Mmax = B->rmap->n; 2180 b->Nmax = B->cmap->n; 2181 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2182 2183 if (!data) { /* petsc-allocated storage */ 2184 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2185 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 2186 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2187 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2188 2189 b->user_alloc = PETSC_FALSE; 2190 } else { /* user-allocated storage */ 2191 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2192 b->v = data; 2193 b->user_alloc = PETSC_TRUE; 2194 } 2195 B->assembled = PETSC_TRUE; 2196 PetscFunctionReturn(0); 2197 } 2198 EXTERN_C_END 2199 2200 #undef __FUNCT__ 2201 #define __FUNCT__ "MatSeqDenseSetLDA" 2202 /*@C 2203 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2204 2205 Input parameter: 2206 + A - the matrix 2207 - lda - the leading dimension 2208 2209 Notes: 2210 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2211 it asserts that the preallocation has a leading dimension (the LDA parameter 2212 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2213 2214 Level: intermediate 2215 2216 .keywords: dense, matrix, LAPACK, BLAS 2217 2218 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2219 2220 @*/ 2221 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2222 { 2223 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2224 2225 PetscFunctionBegin; 2226 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); 2227 b->lda = lda; 2228 b->changelda = PETSC_FALSE; 2229 b->Mmax = PetscMax(b->Mmax,lda); 2230 PetscFunctionReturn(0); 2231 } 2232 2233 /*MC 2234 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2235 2236 Options Database Keys: 2237 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2238 2239 Level: beginner 2240 2241 .seealso: MatCreateSeqDense() 2242 2243 M*/ 2244 2245 EXTERN_C_BEGIN 2246 #undef __FUNCT__ 2247 #define __FUNCT__ "MatCreate_SeqDense" 2248 PetscErrorCode MatCreate_SeqDense(Mat B) 2249 { 2250 Mat_SeqDense *b; 2251 PetscErrorCode ierr; 2252 PetscMPIInt size; 2253 2254 PetscFunctionBegin; 2255 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2256 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2257 2258 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2259 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2260 B->data = (void*)b; 2261 2262 b->pivots = 0; 2263 b->roworiented = PETSC_TRUE; 2264 b->v = 0; 2265 b->changelda = PETSC_FALSE; 2266 2267 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseGetArray_C","MatDenseGetArray_SeqDense",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2268 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatDenseRestoreArray_C","MatDenseRestoreArray_SeqDense",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2269 2270 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_seqaij_C","MatConvert_SeqDense_SeqAIJ",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2271 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2272 "MatGetFactor_seqdense_petsc", 2273 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2274 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2275 "MatSeqDenseSetPreallocation_SeqDense", 2276 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2277 2278 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2279 "MatMatMult_SeqAIJ_SeqDense", 2280 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2281 2282 2283 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2284 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2285 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2287 "MatMatMultNumeric_SeqAIJ_SeqDense", 2288 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2289 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2290 PetscFunctionReturn(0); 2291 } 2292 EXTERN_C_END 2293