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