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 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 PetscInt i, j; 535 PetscBLASInt m, n, _One=1; 536 537 PetscFunctionBegin; 538 m = PetscBLASIntCast(A->rmap->n); 539 n = PetscBLASIntCast(A->cmap->n); 540 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 541 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 542 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 543 #if 0 544 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 545 { 546 PetscScalar *yCheck = new PetscScalar[m]; 547 for(int i = 0; i < m; ++i) { 548 yCheck[i] = 0.0; 549 for(int j = 0; j < n; ++j) { 550 yCheck[i] += v[i*n+j]*x[j]; 551 } 552 if (fabs(y[i] - yCheck[i]) > 1.0e-8) { 553 SETERRQ2(PETSC_ERR_LIB, "DGEMV call is fucked up: %g should be %g", y[i], yCheck[i]); 554 } 555 } 556 delete [] y; 557 } 558 #else 559 for(i = 0; i < m; ++i) { 560 y[i] = 0.0; 561 for(j = 0; j < n; ++j) { 562 y[i] += v[i*n+j]*x[j]; 563 } 564 } 565 #endif 566 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 567 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 568 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 569 PetscFunctionReturn(0); 570 } 571 572 #undef __FUNCT__ 573 #define __FUNCT__ "MatMultAdd_SeqDense" 574 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 575 { 576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 577 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 578 PetscErrorCode ierr; 579 PetscBLASInt m, n, _One=1; 580 581 PetscFunctionBegin; 582 m = PetscBLASIntCast(A->rmap->n); 583 n = PetscBLASIntCast(A->cmap->n); 584 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 585 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 586 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 587 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 588 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 589 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 590 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 591 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 #undef __FUNCT__ 596 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 597 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 598 { 599 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 600 PetscScalar *v = mat->v,*x,*y; 601 PetscErrorCode ierr; 602 PetscBLASInt m, n, _One=1; 603 PetscScalar _DOne=1.0; 604 605 PetscFunctionBegin; 606 m = PetscBLASIntCast(A->rmap->n); 607 n = PetscBLASIntCast(A->cmap->n); 608 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 609 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 610 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 611 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 612 BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 613 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 614 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 615 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 616 PetscFunctionReturn(0); 617 } 618 619 /* -----------------------------------------------------------------*/ 620 #undef __FUNCT__ 621 #define __FUNCT__ "MatGetRow_SeqDense" 622 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 623 { 624 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 625 PetscScalar *v; 626 PetscErrorCode ierr; 627 PetscInt i; 628 629 PetscFunctionBegin; 630 *ncols = A->cmap->n; 631 if (cols) { 632 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);CHKERRQ(ierr); 633 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 634 } 635 if (vals) { 636 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 637 v = mat->v + row; 638 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 639 } 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "MatRestoreRow_SeqDense" 645 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 646 { 647 PetscErrorCode ierr; 648 PetscFunctionBegin; 649 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 650 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 651 PetscFunctionReturn(0); 652 } 653 /* ----------------------------------------------------------------*/ 654 #undef __FUNCT__ 655 #define __FUNCT__ "MatSetValues_SeqDense" 656 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 657 { 658 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 659 PetscInt i,j,idx=0; 660 661 PetscFunctionBegin; 662 if (v) PetscValidScalarPointer(v,6); 663 if (!mat->roworiented) { 664 if (addv == INSERT_VALUES) { 665 for (j=0; j<n; j++) { 666 if (indexn[j] < 0) {idx += m; continue;} 667 #if defined(PETSC_USE_DEBUG) 668 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 669 #endif 670 for (i=0; i<m; i++) { 671 if (indexm[i] < 0) {idx++; continue;} 672 #if defined(PETSC_USE_DEBUG) 673 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 674 #endif 675 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 676 } 677 } 678 } else { 679 for (j=0; j<n; j++) { 680 if (indexn[j] < 0) {idx += m; continue;} 681 #if defined(PETSC_USE_DEBUG) 682 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 683 #endif 684 for (i=0; i<m; i++) { 685 if (indexm[i] < 0) {idx++; continue;} 686 #if defined(PETSC_USE_DEBUG) 687 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 688 #endif 689 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 690 } 691 } 692 } 693 } else { 694 if (addv == INSERT_VALUES) { 695 for (i=0; i<m; i++) { 696 if (indexm[i] < 0) { idx += n; continue;} 697 #if defined(PETSC_USE_DEBUG) 698 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 699 #endif 700 for (j=0; j<n; j++) { 701 if (indexn[j] < 0) { idx++; continue;} 702 #if defined(PETSC_USE_DEBUG) 703 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 704 #endif 705 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 706 } 707 } 708 } else { 709 for (i=0; i<m; i++) { 710 if (indexm[i] < 0) { idx += n; continue;} 711 #if defined(PETSC_USE_DEBUG) 712 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 713 #endif 714 for (j=0; j<n; j++) { 715 if (indexn[j] < 0) { idx++; continue;} 716 #if defined(PETSC_USE_DEBUG) 717 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 718 #endif 719 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 720 } 721 } 722 } 723 } 724 PetscFunctionReturn(0); 725 } 726 727 #undef __FUNCT__ 728 #define __FUNCT__ "MatGetValues_SeqDense" 729 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 730 { 731 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 732 PetscInt i,j; 733 734 PetscFunctionBegin; 735 /* row-oriented output */ 736 for (i=0; i<m; i++) { 737 if (indexm[i] < 0) {v += n;continue;} 738 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 739 for (j=0; j<n; j++) { 740 if (indexn[j] < 0) {v++; continue;} 741 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 742 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 743 } 744 } 745 PetscFunctionReturn(0); 746 } 747 748 /* -----------------------------------------------------------------*/ 749 750 #include "petscsys.h" 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "MatLoad_SeqDense" 754 PetscErrorCode MatLoad_SeqDense(PetscViewer viewer, const MatType type,Mat *A) 755 { 756 Mat_SeqDense *a; 757 Mat B; 758 PetscErrorCode ierr; 759 PetscInt *scols,i,j,nz,header[4]; 760 int fd; 761 PetscMPIInt size; 762 PetscInt *rowlengths = 0,M,N,*cols; 763 PetscScalar *vals,*svals,*v,*w; 764 MPI_Comm comm = ((PetscObject)viewer)->comm; 765 766 PetscFunctionBegin; 767 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 768 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 769 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 770 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 771 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 772 M = header[1]; N = header[2]; nz = header[3]; 773 774 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 775 ierr = MatCreate(comm,A);CHKERRQ(ierr); 776 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 777 ierr = MatSetType(*A,type);CHKERRQ(ierr); 778 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 779 B = *A; 780 a = (Mat_SeqDense*)B->data; 781 v = a->v; 782 /* Allocate some temp space to read in the values and then flip them 783 from row major to column major */ 784 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 785 /* read in nonzero values */ 786 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 787 /* now flip the values and store them in the matrix*/ 788 for (j=0; j<N; j++) { 789 for (i=0; i<M; i++) { 790 *v++ =w[i*N+j]; 791 } 792 } 793 ierr = PetscFree(w);CHKERRQ(ierr); 794 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 795 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 796 } else { 797 /* read row lengths */ 798 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 799 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 800 801 /* create our matrix */ 802 ierr = MatCreate(comm,A);CHKERRQ(ierr); 803 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 804 ierr = MatSetType(*A,type);CHKERRQ(ierr); 805 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 806 B = *A; 807 a = (Mat_SeqDense*)B->data; 808 v = a->v; 809 810 /* read column indices and nonzeros */ 811 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 812 cols = scols; 813 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 814 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 815 vals = svals; 816 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 817 818 /* insert into matrix */ 819 for (i=0; i<M; i++) { 820 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 821 svals += rowlengths[i]; scols += rowlengths[i]; 822 } 823 ierr = PetscFree(vals);CHKERRQ(ierr); 824 ierr = PetscFree(cols);CHKERRQ(ierr); 825 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 826 827 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 828 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 829 } 830 PetscFunctionReturn(0); 831 } 832 833 #include "petscsys.h" 834 835 #undef __FUNCT__ 836 #define __FUNCT__ "MatView_SeqDense_ASCII" 837 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 838 { 839 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 840 PetscErrorCode ierr; 841 PetscInt i,j; 842 const char *name; 843 PetscScalar *v; 844 PetscViewerFormat format; 845 #if defined(PETSC_USE_COMPLEX) 846 PetscTruth allreal = PETSC_TRUE; 847 #endif 848 849 PetscFunctionBegin; 850 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 851 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 852 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 853 PetscFunctionReturn(0); /* do nothing for now */ 854 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 855 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 856 for (i=0; i<A->rmap->n; i++) { 857 v = a->v + i; 858 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 859 for (j=0; j<A->cmap->n; j++) { 860 #if defined(PETSC_USE_COMPLEX) 861 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 862 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 863 } else if (PetscRealPart(*v)) { 864 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 865 } 866 #else 867 if (*v) { 868 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 869 } 870 #endif 871 v += a->lda; 872 } 873 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 874 } 875 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 876 } else { 877 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 878 #if defined(PETSC_USE_COMPLEX) 879 /* determine if matrix has all real values */ 880 v = a->v; 881 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 882 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 883 } 884 #endif 885 if (format == PETSC_VIEWER_ASCII_MATLAB) { 886 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 887 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 888 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 889 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 890 } 891 892 for (i=0; i<A->rmap->n; i++) { 893 v = a->v + i; 894 for (j=0; j<A->cmap->n; j++) { 895 #if defined(PETSC_USE_COMPLEX) 896 if (allreal) { 897 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 898 } else { 899 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 900 } 901 #else 902 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 903 #endif 904 v += a->lda; 905 } 906 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 907 } 908 if (format == PETSC_VIEWER_ASCII_MATLAB) { 909 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 910 } 911 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 912 } 913 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "MatView_SeqDense_Binary" 919 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 920 { 921 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 922 PetscErrorCode ierr; 923 int fd; 924 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 925 PetscScalar *v,*anonz,*vals; 926 PetscViewerFormat format; 927 928 PetscFunctionBegin; 929 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 930 931 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 932 if (format == PETSC_VIEWER_NATIVE) { 933 /* store the matrix as a dense matrix */ 934 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 935 col_lens[0] = MAT_FILE_COOKIE; 936 col_lens[1] = m; 937 col_lens[2] = n; 938 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 939 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 940 ierr = PetscFree(col_lens);CHKERRQ(ierr); 941 942 /* write out matrix, by rows */ 943 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 944 v = a->v; 945 for (j=0; j<n; j++) { 946 for (i=0; i<m; i++) { 947 vals[j + i*n] = *v++; 948 } 949 } 950 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 951 ierr = PetscFree(vals);CHKERRQ(ierr); 952 } else { 953 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 954 col_lens[0] = MAT_FILE_COOKIE; 955 col_lens[1] = m; 956 col_lens[2] = n; 957 col_lens[3] = nz; 958 959 /* store lengths of each row and write (including header) to file */ 960 for (i=0; i<m; i++) col_lens[4+i] = n; 961 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 962 963 /* Possibly should write in smaller increments, not whole matrix at once? */ 964 /* store column indices (zero start index) */ 965 ict = 0; 966 for (i=0; i<m; i++) { 967 for (j=0; j<n; j++) col_lens[ict++] = j; 968 } 969 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 970 ierr = PetscFree(col_lens);CHKERRQ(ierr); 971 972 /* store nonzero values */ 973 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 974 ict = 0; 975 for (i=0; i<m; i++) { 976 v = a->v + i; 977 for (j=0; j<n; j++) { 978 anonz[ict++] = *v; v += a->lda; 979 } 980 } 981 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 982 ierr = PetscFree(anonz);CHKERRQ(ierr); 983 } 984 PetscFunctionReturn(0); 985 } 986 987 #undef __FUNCT__ 988 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 989 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 990 { 991 Mat A = (Mat) Aa; 992 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 993 PetscErrorCode ierr; 994 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 995 PetscScalar *v = a->v; 996 PetscViewer viewer; 997 PetscDraw popup; 998 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 999 PetscViewerFormat format; 1000 1001 PetscFunctionBegin; 1002 1003 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1004 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1005 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1006 1007 /* Loop over matrix elements drawing boxes */ 1008 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1009 /* Blue for negative and Red for positive */ 1010 color = PETSC_DRAW_BLUE; 1011 for(j = 0; j < n; j++) { 1012 x_l = j; 1013 x_r = x_l + 1.0; 1014 for(i = 0; i < m; i++) { 1015 y_l = m - i - 1.0; 1016 y_r = y_l + 1.0; 1017 #if defined(PETSC_USE_COMPLEX) 1018 if (PetscRealPart(v[j*m+i]) > 0.) { 1019 color = PETSC_DRAW_RED; 1020 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1021 color = PETSC_DRAW_BLUE; 1022 } else { 1023 continue; 1024 } 1025 #else 1026 if (v[j*m+i] > 0.) { 1027 color = PETSC_DRAW_RED; 1028 } else if (v[j*m+i] < 0.) { 1029 color = PETSC_DRAW_BLUE; 1030 } else { 1031 continue; 1032 } 1033 #endif 1034 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1035 } 1036 } 1037 } else { 1038 /* use contour shading to indicate magnitude of values */ 1039 /* first determine max of all nonzero values */ 1040 for(i = 0; i < m*n; i++) { 1041 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1042 } 1043 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1044 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1045 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1046 for(j = 0; j < n; j++) { 1047 x_l = j; 1048 x_r = x_l + 1.0; 1049 for(i = 0; i < m; i++) { 1050 y_l = m - i - 1.0; 1051 y_r = y_l + 1.0; 1052 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1053 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1054 } 1055 } 1056 } 1057 PetscFunctionReturn(0); 1058 } 1059 1060 #undef __FUNCT__ 1061 #define __FUNCT__ "MatView_SeqDense_Draw" 1062 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1063 { 1064 PetscDraw draw; 1065 PetscTruth isnull; 1066 PetscReal xr,yr,xl,yl,h,w; 1067 PetscErrorCode ierr; 1068 1069 PetscFunctionBegin; 1070 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1071 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1072 if (isnull) PetscFunctionReturn(0); 1073 1074 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1075 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1076 xr += w; yr += h; xl = -w; yl = -h; 1077 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1078 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1079 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1080 PetscFunctionReturn(0); 1081 } 1082 1083 #undef __FUNCT__ 1084 #define __FUNCT__ "MatView_SeqDense" 1085 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1086 { 1087 PetscErrorCode ierr; 1088 PetscTruth iascii,isbinary,isdraw; 1089 1090 PetscFunctionBegin; 1091 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 1092 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 1093 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 1094 1095 if (iascii) { 1096 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1097 } else if (isbinary) { 1098 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1099 } else if (isdraw) { 1100 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1101 } else { 1102 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1103 } 1104 PetscFunctionReturn(0); 1105 } 1106 1107 #undef __FUNCT__ 1108 #define __FUNCT__ "MatDestroy_SeqDense" 1109 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1110 { 1111 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1112 PetscErrorCode ierr; 1113 1114 PetscFunctionBegin; 1115 #if defined(PETSC_USE_LOG) 1116 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1117 #endif 1118 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1119 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1120 ierr = PetscFree(l);CHKERRQ(ierr); 1121 1122 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1123 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1124 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1125 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1126 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 #undef __FUNCT__ 1131 #define __FUNCT__ "MatTranspose_SeqDense" 1132 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1133 { 1134 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1135 PetscErrorCode ierr; 1136 PetscInt k,j,m,n,M; 1137 PetscScalar *v,tmp; 1138 1139 PetscFunctionBegin; 1140 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1141 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1142 if (m != n) { 1143 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1144 } else { 1145 for (j=0; j<m; j++) { 1146 for (k=0; k<j; k++) { 1147 tmp = v[j + k*M]; 1148 v[j + k*M] = v[k + j*M]; 1149 v[k + j*M] = tmp; 1150 } 1151 } 1152 } 1153 } else { /* out-of-place transpose */ 1154 Mat tmat; 1155 Mat_SeqDense *tmatd; 1156 PetscScalar *v2; 1157 1158 if (reuse == MAT_INITIAL_MATRIX) { 1159 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1160 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1161 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1162 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1163 } else { 1164 tmat = *matout; 1165 } 1166 tmatd = (Mat_SeqDense*)tmat->data; 1167 v = mat->v; v2 = tmatd->v; 1168 for (j=0; j<n; j++) { 1169 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1170 } 1171 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1172 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1173 *matout = tmat; 1174 } 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "MatEqual_SeqDense" 1180 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1181 { 1182 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1183 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1184 PetscInt i,j; 1185 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1186 1187 PetscFunctionBegin; 1188 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1189 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1190 for (i=0; i<A1->rmap->n; i++) { 1191 v1 = mat1->v+i; v2 = mat2->v+i; 1192 for (j=0; j<A1->cmap->n; j++) { 1193 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1194 v1 += mat1->lda; v2 += mat2->lda; 1195 } 1196 } 1197 *flg = PETSC_TRUE; 1198 PetscFunctionReturn(0); 1199 } 1200 1201 #undef __FUNCT__ 1202 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1203 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1204 { 1205 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1206 PetscErrorCode ierr; 1207 PetscInt i,n,len; 1208 PetscScalar *x,zero = 0.0; 1209 1210 PetscFunctionBegin; 1211 ierr = VecSet(v,zero);CHKERRQ(ierr); 1212 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1213 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1214 len = PetscMin(A->rmap->n,A->cmap->n); 1215 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1216 for (i=0; i<len; i++) { 1217 x[i] = mat->v[i*mat->lda + i]; 1218 } 1219 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1220 PetscFunctionReturn(0); 1221 } 1222 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1225 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1226 { 1227 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1228 PetscScalar *l,*r,x,*v; 1229 PetscErrorCode ierr; 1230 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1231 1232 PetscFunctionBegin; 1233 if (ll) { 1234 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1235 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1236 if (m != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1237 for (i=0; i<m; i++) { 1238 x = l[i]; 1239 v = mat->v + i; 1240 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1241 } 1242 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1243 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1244 } 1245 if (rr) { 1246 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1247 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1248 if (n != A->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1249 for (i=0; i<n; i++) { 1250 x = r[i]; 1251 v = mat->v + i*m; 1252 for (j=0; j<m; j++) { (*v++) *= x;} 1253 } 1254 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1255 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1256 } 1257 PetscFunctionReturn(0); 1258 } 1259 1260 #undef __FUNCT__ 1261 #define __FUNCT__ "MatNorm_SeqDense" 1262 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1263 { 1264 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1265 PetscScalar *v = mat->v; 1266 PetscReal sum = 0.0; 1267 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1268 PetscErrorCode ierr; 1269 1270 PetscFunctionBegin; 1271 if (type == NORM_FROBENIUS) { 1272 if (lda>m) { 1273 for (j=0; j<A->cmap->n; j++) { 1274 v = mat->v+j*lda; 1275 for (i=0; i<m; i++) { 1276 #if defined(PETSC_USE_COMPLEX) 1277 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1278 #else 1279 sum += (*v)*(*v); v++; 1280 #endif 1281 } 1282 } 1283 } else { 1284 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1285 #if defined(PETSC_USE_COMPLEX) 1286 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1287 #else 1288 sum += (*v)*(*v); v++; 1289 #endif 1290 } 1291 } 1292 *nrm = sqrt(sum); 1293 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1294 } else if (type == NORM_1) { 1295 *nrm = 0.0; 1296 for (j=0; j<A->cmap->n; j++) { 1297 v = mat->v + j*mat->lda; 1298 sum = 0.0; 1299 for (i=0; i<A->rmap->n; i++) { 1300 sum += PetscAbsScalar(*v); v++; 1301 } 1302 if (sum > *nrm) *nrm = sum; 1303 } 1304 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1305 } else if (type == NORM_INFINITY) { 1306 *nrm = 0.0; 1307 for (j=0; j<A->rmap->n; j++) { 1308 v = mat->v + j; 1309 sum = 0.0; 1310 for (i=0; i<A->cmap->n; i++) { 1311 sum += PetscAbsScalar(*v); v += mat->lda; 1312 } 1313 if (sum > *nrm) *nrm = sum; 1314 } 1315 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1316 } else { 1317 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1318 } 1319 PetscFunctionReturn(0); 1320 } 1321 1322 #undef __FUNCT__ 1323 #define __FUNCT__ "MatSetOption_SeqDense" 1324 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1325 { 1326 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1327 PetscErrorCode ierr; 1328 1329 PetscFunctionBegin; 1330 switch (op) { 1331 case MAT_ROW_ORIENTED: 1332 aij->roworiented = flg; 1333 break; 1334 case MAT_NEW_NONZERO_LOCATIONS: 1335 case MAT_NEW_NONZERO_LOCATION_ERR: 1336 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1337 case MAT_NEW_DIAGONALS: 1338 case MAT_IGNORE_OFF_PROC_ENTRIES: 1339 case MAT_USE_HASH_TABLE: 1340 case MAT_SYMMETRIC: 1341 case MAT_STRUCTURALLY_SYMMETRIC: 1342 case MAT_HERMITIAN: 1343 case MAT_SYMMETRY_ETERNAL: 1344 case MAT_IGNORE_LOWER_TRIANGULAR: 1345 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1346 break; 1347 default: 1348 SETERRQ1(PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1349 } 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #undef __FUNCT__ 1354 #define __FUNCT__ "MatZeroEntries_SeqDense" 1355 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1356 { 1357 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1358 PetscErrorCode ierr; 1359 PetscInt lda=l->lda,m=A->rmap->n,j; 1360 1361 PetscFunctionBegin; 1362 if (lda>m) { 1363 for (j=0; j<A->cmap->n; j++) { 1364 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1365 } 1366 } else { 1367 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1368 } 1369 PetscFunctionReturn(0); 1370 } 1371 1372 #undef __FUNCT__ 1373 #define __FUNCT__ "MatZeroRows_SeqDense" 1374 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1375 { 1376 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1377 PetscInt n = A->cmap->n,i,j; 1378 PetscScalar *slot; 1379 1380 PetscFunctionBegin; 1381 for (i=0; i<N; i++) { 1382 slot = l->v + rows[i]; 1383 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1384 } 1385 if (diag != 0.0) { 1386 for (i=0; i<N; i++) { 1387 slot = l->v + (n+1)*rows[i]; 1388 *slot = diag; 1389 } 1390 } 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "MatGetArray_SeqDense" 1396 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1397 { 1398 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1399 1400 PetscFunctionBegin; 1401 if (mat->lda != A->rmap->n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1402 *array = mat->v; 1403 PetscFunctionReturn(0); 1404 } 1405 1406 #undef __FUNCT__ 1407 #define __FUNCT__ "MatRestoreArray_SeqDense" 1408 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1409 { 1410 PetscFunctionBegin; 1411 *array = 0; /* user cannot accidently use the array later */ 1412 PetscFunctionReturn(0); 1413 } 1414 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1417 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1418 { 1419 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1420 PetscErrorCode ierr; 1421 PetscInt i,j,nrows,ncols; 1422 const PetscInt *irow,*icol; 1423 PetscScalar *av,*bv,*v = mat->v; 1424 Mat newmat; 1425 1426 PetscFunctionBegin; 1427 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1428 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1429 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1430 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1431 1432 /* Check submatrixcall */ 1433 if (scall == MAT_REUSE_MATRIX) { 1434 PetscInt n_cols,n_rows; 1435 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1436 if (n_rows != nrows || n_cols != ncols) { 1437 /* resize the result result matrix to match number of requested rows/columns */ 1438 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1439 } 1440 newmat = *B; 1441 } else { 1442 /* Create and fill new matrix */ 1443 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1444 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1445 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1446 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1447 } 1448 1449 /* Now extract the data pointers and do the copy,column at a time */ 1450 bv = ((Mat_SeqDense*)newmat->data)->v; 1451 1452 for (i=0; i<ncols; i++) { 1453 av = v + mat->lda*icol[i]; 1454 for (j=0; j<nrows; j++) { 1455 *bv++ = av[irow[j]]; 1456 } 1457 } 1458 1459 /* Assemble the matrices so that the correct flags are set */ 1460 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1461 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1462 1463 /* Free work space */ 1464 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1465 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1466 *B = newmat; 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1472 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1473 { 1474 PetscErrorCode ierr; 1475 PetscInt i; 1476 1477 PetscFunctionBegin; 1478 if (scall == MAT_INITIAL_MATRIX) { 1479 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1480 } 1481 1482 for (i=0; i<n; i++) { 1483 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1484 } 1485 PetscFunctionReturn(0); 1486 } 1487 1488 #undef __FUNCT__ 1489 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1490 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1491 { 1492 PetscFunctionBegin; 1493 PetscFunctionReturn(0); 1494 } 1495 1496 #undef __FUNCT__ 1497 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1498 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1499 { 1500 PetscFunctionBegin; 1501 PetscFunctionReturn(0); 1502 } 1503 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "MatCopy_SeqDense" 1506 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1507 { 1508 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1509 PetscErrorCode ierr; 1510 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1511 1512 PetscFunctionBegin; 1513 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1514 if (A->ops->copy != B->ops->copy) { 1515 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1516 PetscFunctionReturn(0); 1517 } 1518 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1519 if (lda1>m || lda2>m) { 1520 for (j=0; j<n; j++) { 1521 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1522 } 1523 } else { 1524 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1525 } 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #undef __FUNCT__ 1530 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1531 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1532 { 1533 PetscErrorCode ierr; 1534 1535 PetscFunctionBegin; 1536 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatSetSizes_SeqDense" 1542 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1543 { 1544 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1545 PetscFunctionBegin; 1546 /* this will not be called before lda, Mmax, and Nmax have been set */ 1547 m = PetscMax(m,M); 1548 n = PetscMax(n,N); 1549 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); 1550 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); 1551 A->rmap->n = A->rmap->n = m; 1552 A->cmap->n = A->cmap->N = n; 1553 if (a->changelda) a->lda = m; 1554 PetscFunctionReturn(0); 1555 } 1556 1557 1558 /* ----------------------------------------------------------------*/ 1559 #undef __FUNCT__ 1560 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1561 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1562 { 1563 PetscErrorCode ierr; 1564 1565 PetscFunctionBegin; 1566 if (scall == MAT_INITIAL_MATRIX){ 1567 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1568 } 1569 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 #undef __FUNCT__ 1574 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1575 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1576 { 1577 PetscErrorCode ierr; 1578 PetscInt m=A->rmap->n,n=B->cmap->n; 1579 Mat Cmat; 1580 1581 PetscFunctionBegin; 1582 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); 1583 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1584 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1585 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1586 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1587 Cmat->assembled = PETSC_TRUE; 1588 *C = Cmat; 1589 PetscFunctionReturn(0); 1590 } 1591 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1594 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1595 { 1596 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1597 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1598 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1599 PetscBLASInt m,n,k; 1600 PetscScalar _DOne=1.0,_DZero=0.0; 1601 1602 PetscFunctionBegin; 1603 m = PetscBLASIntCast(A->rmap->n); 1604 n = PetscBLASIntCast(B->cmap->n); 1605 k = PetscBLASIntCast(A->cmap->n); 1606 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1607 PetscFunctionReturn(0); 1608 } 1609 1610 #undef __FUNCT__ 1611 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1612 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1613 { 1614 PetscErrorCode ierr; 1615 1616 PetscFunctionBegin; 1617 if (scall == MAT_INITIAL_MATRIX){ 1618 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1619 } 1620 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1621 PetscFunctionReturn(0); 1622 } 1623 1624 #undef __FUNCT__ 1625 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1626 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1627 { 1628 PetscErrorCode ierr; 1629 PetscInt m=A->cmap->n,n=B->cmap->n; 1630 Mat Cmat; 1631 1632 PetscFunctionBegin; 1633 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); 1634 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1635 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1636 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1637 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1638 Cmat->assembled = PETSC_TRUE; 1639 *C = Cmat; 1640 PetscFunctionReturn(0); 1641 } 1642 1643 #undef __FUNCT__ 1644 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1645 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1646 { 1647 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1648 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1649 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1650 PetscBLASInt m,n,k; 1651 PetscScalar _DOne=1.0,_DZero=0.0; 1652 1653 PetscFunctionBegin; 1654 m = PetscBLASIntCast(A->cmap->n); 1655 n = PetscBLASIntCast(B->cmap->n); 1656 k = PetscBLASIntCast(A->rmap->n); 1657 /* 1658 Note the m and n arguments below are the number rows and columns of A', not A! 1659 */ 1660 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1661 PetscFunctionReturn(0); 1662 } 1663 1664 #undef __FUNCT__ 1665 #define __FUNCT__ "MatGetRowMax_SeqDense" 1666 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1667 { 1668 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1669 PetscErrorCode ierr; 1670 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1671 PetscScalar *x; 1672 MatScalar *aa = a->v; 1673 1674 PetscFunctionBegin; 1675 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1676 1677 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1678 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1679 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1680 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1681 for (i=0; i<m; i++) { 1682 x[i] = aa[i]; if (idx) idx[i] = 0; 1683 for (j=1; j<n; j++){ 1684 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1685 } 1686 } 1687 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1688 PetscFunctionReturn(0); 1689 } 1690 1691 #undef __FUNCT__ 1692 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1693 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1694 { 1695 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1696 PetscErrorCode ierr; 1697 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1698 PetscScalar *x; 1699 PetscReal atmp; 1700 MatScalar *aa = a->v; 1701 1702 PetscFunctionBegin; 1703 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1704 1705 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1706 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1707 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1708 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1709 for (i=0; i<m; i++) { 1710 x[i] = PetscAbsScalar(aa[i]); 1711 for (j=1; j<n; j++){ 1712 atmp = PetscAbsScalar(aa[i+m*j]); 1713 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1714 } 1715 } 1716 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1717 PetscFunctionReturn(0); 1718 } 1719 1720 #undef __FUNCT__ 1721 #define __FUNCT__ "MatGetRowMin_SeqDense" 1722 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1723 { 1724 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1725 PetscErrorCode ierr; 1726 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1727 PetscScalar *x; 1728 MatScalar *aa = a->v; 1729 1730 PetscFunctionBegin; 1731 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1732 1733 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1734 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1735 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1736 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1737 for (i=0; i<m; i++) { 1738 x[i] = aa[i]; if (idx) idx[i] = 0; 1739 for (j=1; j<n; j++){ 1740 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1741 } 1742 } 1743 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1744 PetscFunctionReturn(0); 1745 } 1746 1747 #undef __FUNCT__ 1748 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1749 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1750 { 1751 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1752 PetscErrorCode ierr; 1753 PetscScalar *x; 1754 1755 PetscFunctionBegin; 1756 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1757 1758 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1759 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1760 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1761 PetscFunctionReturn(0); 1762 } 1763 1764 /* -------------------------------------------------------------------*/ 1765 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1766 MatGetRow_SeqDense, 1767 MatRestoreRow_SeqDense, 1768 MatMult_SeqDense, 1769 /* 4*/ MatMultAdd_SeqDense, 1770 MatMultTranspose_SeqDense, 1771 MatMultTransposeAdd_SeqDense, 1772 0, 1773 0, 1774 0, 1775 /*10*/ 0, 1776 MatLUFactor_SeqDense, 1777 MatCholeskyFactor_SeqDense, 1778 MatSOR_SeqDense, 1779 MatTranspose_SeqDense, 1780 /*15*/ MatGetInfo_SeqDense, 1781 MatEqual_SeqDense, 1782 MatGetDiagonal_SeqDense, 1783 MatDiagonalScale_SeqDense, 1784 MatNorm_SeqDense, 1785 /*20*/ MatAssemblyBegin_SeqDense, 1786 MatAssemblyEnd_SeqDense, 1787 MatSetOption_SeqDense, 1788 MatZeroEntries_SeqDense, 1789 /*24*/ MatZeroRows_SeqDense, 1790 0, 1791 0, 1792 0, 1793 0, 1794 /*29*/ MatSetUpPreallocation_SeqDense, 1795 0, 1796 0, 1797 MatGetArray_SeqDense, 1798 MatRestoreArray_SeqDense, 1799 /*34*/ MatDuplicate_SeqDense, 1800 0, 1801 0, 1802 0, 1803 0, 1804 /*39*/ MatAXPY_SeqDense, 1805 MatGetSubMatrices_SeqDense, 1806 0, 1807 MatGetValues_SeqDense, 1808 MatCopy_SeqDense, 1809 /*44*/ MatGetRowMax_SeqDense, 1810 MatScale_SeqDense, 1811 0, 1812 0, 1813 0, 1814 /*49*/ 0, 1815 0, 1816 0, 1817 0, 1818 0, 1819 /*54*/ 0, 1820 0, 1821 0, 1822 0, 1823 0, 1824 /*59*/ 0, 1825 MatDestroy_SeqDense, 1826 MatView_SeqDense, 1827 0, 1828 0, 1829 /*64*/ 0, 1830 0, 1831 0, 1832 0, 1833 0, 1834 /*69*/ MatGetRowMaxAbs_SeqDense, 1835 0, 1836 0, 1837 0, 1838 0, 1839 /*74*/ 0, 1840 0, 1841 0, 1842 0, 1843 0, 1844 /*79*/ 0, 1845 0, 1846 0, 1847 0, 1848 /*83*/ MatLoad_SeqDense, 1849 0, 1850 MatIsHermitian_SeqDense, 1851 0, 1852 0, 1853 0, 1854 /*89*/ MatMatMult_SeqDense_SeqDense, 1855 MatMatMultSymbolic_SeqDense_SeqDense, 1856 MatMatMultNumeric_SeqDense_SeqDense, 1857 0, 1858 0, 1859 /*94*/ 0, 1860 MatMatMultTranspose_SeqDense_SeqDense, 1861 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1862 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1863 0, 1864 /*99*/ 0, 1865 0, 1866 0, 1867 0, 1868 MatSetSizes_SeqDense, 1869 0, 1870 0, 1871 0, 1872 0, 1873 0, 1874 /*109*/0, 1875 0, 1876 MatGetRowMin_SeqDense, 1877 MatGetColumnVector_SeqDense 1878 }; 1879 1880 #undef __FUNCT__ 1881 #define __FUNCT__ "MatCreateSeqDense" 1882 /*@C 1883 MatCreateSeqDense - Creates a sequential dense matrix that 1884 is stored in column major order (the usual Fortran 77 manner). Many 1885 of the matrix operations use the BLAS and LAPACK routines. 1886 1887 Collective on MPI_Comm 1888 1889 Input Parameters: 1890 + comm - MPI communicator, set to PETSC_COMM_SELF 1891 . m - number of rows 1892 . n - number of columns 1893 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1894 to control all matrix memory allocation. 1895 1896 Output Parameter: 1897 . A - the matrix 1898 1899 Notes: 1900 The data input variable is intended primarily for Fortran programmers 1901 who wish to allocate their own matrix memory space. Most users should 1902 set data=PETSC_NULL. 1903 1904 Level: intermediate 1905 1906 .keywords: dense, matrix, LAPACK, BLAS 1907 1908 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1909 @*/ 1910 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1911 { 1912 PetscErrorCode ierr; 1913 1914 PetscFunctionBegin; 1915 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1916 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1917 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1918 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1919 PetscFunctionReturn(0); 1920 } 1921 1922 #undef __FUNCT__ 1923 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1924 /*@C 1925 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1926 1927 Collective on MPI_Comm 1928 1929 Input Parameters: 1930 + A - the matrix 1931 - data - the array (or PETSC_NULL) 1932 1933 Notes: 1934 The data input variable is intended primarily for Fortran programmers 1935 who wish to allocate their own matrix memory space. Most users should 1936 need not call this routine. 1937 1938 Level: intermediate 1939 1940 .keywords: dense, matrix, LAPACK, BLAS 1941 1942 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues(), MatSeqDenseSetLDA() 1943 1944 @*/ 1945 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1946 { 1947 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1948 1949 PetscFunctionBegin; 1950 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1951 if (f) { 1952 ierr = (*f)(B,data);CHKERRQ(ierr); 1953 } 1954 PetscFunctionReturn(0); 1955 } 1956 1957 EXTERN_C_BEGIN 1958 #undef __FUNCT__ 1959 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1960 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1961 { 1962 Mat_SeqDense *b; 1963 PetscErrorCode ierr; 1964 1965 PetscFunctionBegin; 1966 B->preallocated = PETSC_TRUE; 1967 b = (Mat_SeqDense*)B->data; 1968 if (b->lda <= 0) b->lda = B->rmap->n; 1969 if (!data) { /* petsc-allocated storage */ 1970 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1971 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1972 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1973 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1974 b->user_alloc = PETSC_FALSE; 1975 } else { /* user-allocated storage */ 1976 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1977 b->v = data; 1978 b->user_alloc = PETSC_TRUE; 1979 } 1980 B->assembled = PETSC_TRUE; 1981 PetscFunctionReturn(0); 1982 } 1983 EXTERN_C_END 1984 1985 #undef __FUNCT__ 1986 #define __FUNCT__ "MatSeqDenseSetLDA" 1987 /*@C 1988 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1989 1990 Input parameter: 1991 + A - the matrix 1992 - lda - the leading dimension 1993 1994 Notes: 1995 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 1996 it asserts that the preallocation has a leading dimension (the LDA parameter 1997 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1998 1999 Level: intermediate 2000 2001 .keywords: dense, matrix, LAPACK, BLAS 2002 2003 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2004 2005 @*/ 2006 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 2007 { 2008 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2009 2010 PetscFunctionBegin; 2011 if (lda < B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2012 b->lda = lda; 2013 b->changelda = PETSC_FALSE; 2014 b->Mmax = PetscMax(b->Mmax,lda); 2015 PetscFunctionReturn(0); 2016 } 2017 2018 /*MC 2019 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2020 2021 Options Database Keys: 2022 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2023 2024 Level: beginner 2025 2026 .seealso: MatCreateSeqDense() 2027 2028 M*/ 2029 2030 EXTERN_C_BEGIN 2031 #undef __FUNCT__ 2032 #define __FUNCT__ "MatCreate_SeqDense" 2033 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2034 { 2035 Mat_SeqDense *b; 2036 PetscErrorCode ierr; 2037 PetscMPIInt size; 2038 2039 PetscFunctionBegin; 2040 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2041 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2042 2043 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2044 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2045 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2046 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2047 2048 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2049 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2050 B->mapping = 0; 2051 B->data = (void*)b; 2052 2053 2054 b->pivots = 0; 2055 b->roworiented = PETSC_TRUE; 2056 b->v = 0; 2057 b->lda = B->rmap->n; 2058 b->changelda = PETSC_FALSE; 2059 b->Mmax = B->rmap->n; 2060 b->Nmax = B->cmap->n; 2061 2062 2063 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2064 "MatGetFactor_seqdense_petsc", 2065 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2066 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2067 "MatSeqDenseSetPreallocation_SeqDense", 2068 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2069 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2070 "MatMatMult_SeqAIJ_SeqDense", 2071 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2072 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2073 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2074 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2075 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2076 "MatMatMultNumeric_SeqAIJ_SeqDense", 2077 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2078 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2079 PetscFunctionReturn(0); 2080 } 2081 2082 2083 EXTERN_C_END 2084