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