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