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