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