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