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