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,nrows,ncols; 1413 const PetscInt *irow,*icol; 1414 PetscScalar *av,*bv,*v = mat->v; 1415 Mat newmat; 1416 1417 PetscFunctionBegin; 1418 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1419 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1420 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1421 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1422 1423 /* Check submatrixcall */ 1424 if (scall == MAT_REUSE_MATRIX) { 1425 PetscInt n_cols,n_rows; 1426 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1427 if (n_rows != nrows || n_cols != ncols) { 1428 /* resize the result result matrix to match number of requested rows/columns */ 1429 ierr = MatSetSizes(*B,nrows,nrows,nrows,nrows);CHKERRQ(ierr); 1430 } 1431 newmat = *B; 1432 } else { 1433 /* Create and fill new matrix */ 1434 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1435 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1436 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1437 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1438 } 1439 1440 /* Now extract the data pointers and do the copy,column at a time */ 1441 bv = ((Mat_SeqDense*)newmat->data)->v; 1442 1443 for (i=0; i<ncols; i++) { 1444 av = v + mat->lda*icol[i]; 1445 for (j=0; j<nrows; j++) { 1446 *bv++ = av[irow[j]]; 1447 } 1448 } 1449 1450 /* Assemble the matrices so that the correct flags are set */ 1451 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1452 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1453 1454 /* Free work space */ 1455 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1456 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1457 *B = newmat; 1458 PetscFunctionReturn(0); 1459 } 1460 1461 #undef __FUNCT__ 1462 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1463 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1464 { 1465 PetscErrorCode ierr; 1466 PetscInt i; 1467 1468 PetscFunctionBegin; 1469 if (scall == MAT_INITIAL_MATRIX) { 1470 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1471 } 1472 1473 for (i=0; i<n; i++) { 1474 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1475 } 1476 PetscFunctionReturn(0); 1477 } 1478 1479 #undef __FUNCT__ 1480 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1481 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1482 { 1483 PetscFunctionBegin; 1484 PetscFunctionReturn(0); 1485 } 1486 1487 #undef __FUNCT__ 1488 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1489 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1490 { 1491 PetscFunctionBegin; 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatCopy_SeqDense" 1497 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1498 { 1499 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1500 PetscErrorCode ierr; 1501 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1502 1503 PetscFunctionBegin; 1504 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1505 if (A->ops->copy != B->ops->copy) { 1506 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1507 PetscFunctionReturn(0); 1508 } 1509 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1510 if (lda1>m || lda2>m) { 1511 for (j=0; j<n; j++) { 1512 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1513 } 1514 } else { 1515 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1516 } 1517 PetscFunctionReturn(0); 1518 } 1519 1520 #undef __FUNCT__ 1521 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1522 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1523 { 1524 PetscErrorCode ierr; 1525 1526 PetscFunctionBegin; 1527 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1528 PetscFunctionReturn(0); 1529 } 1530 1531 #undef __FUNCT__ 1532 #define __FUNCT__ "MatSetSizes_SeqDense" 1533 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1534 { 1535 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1536 PetscErrorCode ierr; 1537 PetscFunctionBegin; 1538 /* this will not be called before lda, Mmax, and Nmax have been set */ 1539 m = PetscMax(m,M); 1540 n = PetscMax(n,N); 1541 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); 1542 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); 1543 A->rmap->n = A->rmap->n = m; 1544 A->cmap->n = A->cmap->N = n; 1545 if (a->changelda) a->lda = m; 1546 ierr = PetscMemzero(a->v,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 1547 PetscFunctionReturn(0); 1548 } 1549 1550 1551 /* ----------------------------------------------------------------*/ 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1554 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1555 { 1556 PetscErrorCode ierr; 1557 1558 PetscFunctionBegin; 1559 if (scall == MAT_INITIAL_MATRIX){ 1560 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1561 } 1562 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1563 PetscFunctionReturn(0); 1564 } 1565 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1568 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1569 { 1570 PetscErrorCode ierr; 1571 PetscInt m=A->rmap->n,n=B->cmap->n; 1572 Mat Cmat; 1573 1574 PetscFunctionBegin; 1575 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); 1576 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1577 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1578 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1579 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1580 Cmat->assembled = PETSC_TRUE; 1581 *C = Cmat; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 #undef __FUNCT__ 1586 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1587 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1588 { 1589 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1590 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1591 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1592 PetscBLASInt m,n,k; 1593 PetscScalar _DOne=1.0,_DZero=0.0; 1594 1595 PetscFunctionBegin; 1596 m = PetscBLASIntCast(A->rmap->n); 1597 n = PetscBLASIntCast(B->cmap->n); 1598 k = PetscBLASIntCast(A->cmap->n); 1599 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1600 PetscFunctionReturn(0); 1601 } 1602 1603 #undef __FUNCT__ 1604 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1605 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1606 { 1607 PetscErrorCode ierr; 1608 1609 PetscFunctionBegin; 1610 if (scall == MAT_INITIAL_MATRIX){ 1611 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1612 } 1613 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1619 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1620 { 1621 PetscErrorCode ierr; 1622 PetscInt m=A->cmap->n,n=B->cmap->n; 1623 Mat Cmat; 1624 1625 PetscFunctionBegin; 1626 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); 1627 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1628 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1629 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1630 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1631 Cmat->assembled = PETSC_TRUE; 1632 *C = Cmat; 1633 PetscFunctionReturn(0); 1634 } 1635 1636 #undef __FUNCT__ 1637 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1638 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1639 { 1640 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1641 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1642 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1643 PetscBLASInt m,n,k; 1644 PetscScalar _DOne=1.0,_DZero=0.0; 1645 1646 PetscFunctionBegin; 1647 m = PetscBLASIntCast(A->cmap->n); 1648 n = PetscBLASIntCast(B->cmap->n); 1649 k = PetscBLASIntCast(A->rmap->n); 1650 /* 1651 Note the m and n arguments below are the number rows and columns of A', not A! 1652 */ 1653 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "MatGetRowMax_SeqDense" 1659 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1660 { 1661 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1662 PetscErrorCode ierr; 1663 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1664 PetscScalar *x; 1665 MatScalar *aa = a->v; 1666 1667 PetscFunctionBegin; 1668 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1669 1670 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1671 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1672 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1673 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1674 for (i=0; i<m; i++) { 1675 x[i] = aa[i]; if (idx) idx[i] = 0; 1676 for (j=1; j<n; j++){ 1677 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1678 } 1679 } 1680 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1686 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1687 { 1688 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1689 PetscErrorCode ierr; 1690 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1691 PetscScalar *x; 1692 PetscReal atmp; 1693 MatScalar *aa = a->v; 1694 1695 PetscFunctionBegin; 1696 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1697 1698 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1699 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1700 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1701 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1702 for (i=0; i<m; i++) { 1703 x[i] = PetscAbsScalar(aa[i]); 1704 for (j=1; j<n; j++){ 1705 atmp = PetscAbsScalar(aa[i+m*j]); 1706 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1707 } 1708 } 1709 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1710 PetscFunctionReturn(0); 1711 } 1712 1713 #undef __FUNCT__ 1714 #define __FUNCT__ "MatGetRowMin_SeqDense" 1715 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1716 { 1717 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1718 PetscErrorCode ierr; 1719 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1720 PetscScalar *x; 1721 MatScalar *aa = a->v; 1722 1723 PetscFunctionBegin; 1724 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1725 1726 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1727 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1728 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1729 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1730 for (i=0; i<m; i++) { 1731 x[i] = aa[i]; if (idx) idx[i] = 0; 1732 for (j=1; j<n; j++){ 1733 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1734 } 1735 } 1736 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1737 PetscFunctionReturn(0); 1738 } 1739 1740 #undef __FUNCT__ 1741 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1742 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1743 { 1744 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1745 PetscErrorCode ierr; 1746 PetscScalar *x; 1747 1748 PetscFunctionBegin; 1749 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1750 1751 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1752 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1753 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1754 PetscFunctionReturn(0); 1755 } 1756 1757 /* -------------------------------------------------------------------*/ 1758 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1759 MatGetRow_SeqDense, 1760 MatRestoreRow_SeqDense, 1761 MatMult_SeqDense, 1762 /* 4*/ MatMultAdd_SeqDense, 1763 MatMultTranspose_SeqDense, 1764 MatMultTransposeAdd_SeqDense, 1765 0, 1766 0, 1767 0, 1768 /*10*/ 0, 1769 MatLUFactor_SeqDense, 1770 MatCholeskyFactor_SeqDense, 1771 MatRelax_SeqDense, 1772 MatTranspose_SeqDense, 1773 /*15*/ MatGetInfo_SeqDense, 1774 MatEqual_SeqDense, 1775 MatGetDiagonal_SeqDense, 1776 MatDiagonalScale_SeqDense, 1777 MatNorm_SeqDense, 1778 /*20*/ MatAssemblyBegin_SeqDense, 1779 MatAssemblyEnd_SeqDense, 1780 0, 1781 MatSetOption_SeqDense, 1782 MatZeroEntries_SeqDense, 1783 /*25*/ MatZeroRows_SeqDense, 1784 0, 1785 0, 1786 0, 1787 0, 1788 /*30*/ MatSetUpPreallocation_SeqDense, 1789 0, 1790 0, 1791 MatGetArray_SeqDense, 1792 MatRestoreArray_SeqDense, 1793 /*35*/ MatDuplicate_SeqDense, 1794 0, 1795 0, 1796 0, 1797 0, 1798 /*40*/ MatAXPY_SeqDense, 1799 MatGetSubMatrices_SeqDense, 1800 0, 1801 MatGetValues_SeqDense, 1802 MatCopy_SeqDense, 1803 /*45*/ MatGetRowMax_SeqDense, 1804 MatScale_SeqDense, 1805 0, 1806 0, 1807 0, 1808 /*50*/ 0, 1809 0, 1810 0, 1811 0, 1812 0, 1813 /*55*/ 0, 1814 0, 1815 0, 1816 0, 1817 0, 1818 /*60*/ 0, 1819 MatDestroy_SeqDense, 1820 MatView_SeqDense, 1821 0, 1822 0, 1823 /*65*/ 0, 1824 0, 1825 0, 1826 0, 1827 0, 1828 /*70*/ MatGetRowMaxAbs_SeqDense, 1829 0, 1830 0, 1831 0, 1832 0, 1833 /*75*/ 0, 1834 0, 1835 0, 1836 0, 1837 0, 1838 /*80*/ 0, 1839 0, 1840 0, 1841 0, 1842 /*84*/ MatLoad_SeqDense, 1843 0, 1844 MatIsHermitian_SeqDense, 1845 0, 1846 0, 1847 0, 1848 /*90*/ MatMatMult_SeqDense_SeqDense, 1849 MatMatMultSymbolic_SeqDense_SeqDense, 1850 MatMatMultNumeric_SeqDense_SeqDense, 1851 0, 1852 0, 1853 /*95*/ 0, 1854 MatMatMultTranspose_SeqDense_SeqDense, 1855 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1856 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1857 0, 1858 /*100*/0, 1859 0, 1860 0, 1861 0, 1862 MatSetSizes_SeqDense, 1863 0, 1864 0, 1865 0, 1866 0, 1867 0, 1868 /*110*/0, 1869 0, 1870 MatGetRowMin_SeqDense, 1871 MatGetColumnVector_SeqDense 1872 }; 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatCreateSeqDense" 1876 /*@C 1877 MatCreateSeqDense - Creates a sequential dense matrix that 1878 is stored in column major order (the usual Fortran 77 manner). Many 1879 of the matrix operations use the BLAS and LAPACK routines. 1880 1881 Collective on MPI_Comm 1882 1883 Input Parameters: 1884 + comm - MPI communicator, set to PETSC_COMM_SELF 1885 . m - number of rows 1886 . n - number of columns 1887 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1888 to control all matrix memory allocation. 1889 1890 Output Parameter: 1891 . A - the matrix 1892 1893 Notes: 1894 The data input variable is intended primarily for Fortran programmers 1895 who wish to allocate their own matrix memory space. Most users should 1896 set data=PETSC_NULL. 1897 1898 Level: intermediate 1899 1900 .keywords: dense, matrix, LAPACK, BLAS 1901 1902 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1903 @*/ 1904 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1905 { 1906 PetscErrorCode ierr; 1907 1908 PetscFunctionBegin; 1909 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1910 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1911 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1912 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1913 PetscFunctionReturn(0); 1914 } 1915 1916 #undef __FUNCT__ 1917 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1918 /*@C 1919 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1920 1921 Collective on MPI_Comm 1922 1923 Input Parameters: 1924 + A - the matrix 1925 - data - the array (or PETSC_NULL) 1926 1927 Notes: 1928 The data input variable is intended primarily for Fortran programmers 1929 who wish to allocate their own matrix memory space. Most users should 1930 need not call this routine. 1931 1932 Level: intermediate 1933 1934 .keywords: dense, matrix, LAPACK, BLAS 1935 1936 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1937 @*/ 1938 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1939 { 1940 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1941 1942 PetscFunctionBegin; 1943 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1944 if (f) { 1945 ierr = (*f)(B,data);CHKERRQ(ierr); 1946 } 1947 PetscFunctionReturn(0); 1948 } 1949 1950 EXTERN_C_BEGIN 1951 #undef __FUNCT__ 1952 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1953 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1954 { 1955 Mat_SeqDense *b; 1956 PetscErrorCode ierr; 1957 1958 PetscFunctionBegin; 1959 B->preallocated = PETSC_TRUE; 1960 b = (Mat_SeqDense*)B->data; 1961 if (b->lda <= 0) b->lda = B->rmap->n; 1962 if (!data) { /* petsc-allocated storage */ 1963 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1964 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1965 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1966 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1967 b->user_alloc = PETSC_FALSE; 1968 } else { /* user-allocated storage */ 1969 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1970 b->v = data; 1971 b->user_alloc = PETSC_TRUE; 1972 } 1973 PetscFunctionReturn(0); 1974 } 1975 EXTERN_C_END 1976 1977 #undef __FUNCT__ 1978 #define __FUNCT__ "MatSeqDenseSetLDA" 1979 /*@C 1980 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1981 1982 Input parameter: 1983 + A - the matrix 1984 - lda - the leading dimension 1985 1986 Notes: 1987 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1988 it asserts that the preallocation has a leading dimension (the LDA parameter 1989 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1990 1991 Level: intermediate 1992 1993 .keywords: dense, matrix, LAPACK, BLAS 1994 1995 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 1996 @*/ 1997 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 1998 { 1999 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2000 2001 PetscFunctionBegin; 2002 if (lda < B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2003 b->lda = lda; 2004 b->changelda = PETSC_FALSE; 2005 b->Mmax = PetscMax(b->Mmax,lda); 2006 PetscFunctionReturn(0); 2007 } 2008 2009 /*MC 2010 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2011 2012 Options Database Keys: 2013 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2014 2015 Level: beginner 2016 2017 .seealso: MatCreateSeqDense() 2018 2019 M*/ 2020 2021 EXTERN_C_BEGIN 2022 #undef __FUNCT__ 2023 #define __FUNCT__ "MatCreate_SeqDense" 2024 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2025 { 2026 Mat_SeqDense *b; 2027 PetscErrorCode ierr; 2028 PetscMPIInt size; 2029 2030 PetscFunctionBegin; 2031 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2032 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2033 2034 B->rmap->bs = B->cmap->bs = 1; 2035 ierr = PetscMapSetUp(B->rmap);CHKERRQ(ierr); 2036 ierr = PetscMapSetUp(B->cmap);CHKERRQ(ierr); 2037 2038 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2039 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2040 B->mapping = 0; 2041 B->data = (void*)b; 2042 2043 2044 b->pivots = 0; 2045 b->roworiented = PETSC_TRUE; 2046 b->v = 0; 2047 b->lda = B->rmap->n; 2048 b->changelda = PETSC_FALSE; 2049 b->Mmax = B->rmap->n; 2050 b->Nmax = B->cmap->n; 2051 2052 2053 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqdense_petsc_C", 2054 "MatGetFactor_seqdense_petsc", 2055 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2056 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2057 "MatSeqDenseSetPreallocation_SeqDense", 2058 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2059 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2060 "MatMatMult_SeqAIJ_SeqDense", 2061 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2062 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2063 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2064 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2065 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2066 "MatMatMultNumeric_SeqAIJ_SeqDense", 2067 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2068 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2069 PetscFunctionReturn(0); 2070 } 2071 2072 2073 EXTERN_C_END 2074