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