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