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,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) { /* 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 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1106 ierr = MatSetSizes(tmat,A->cmap.n,A->rmap.n,A->cmap.n,A->rmap.n);CHKERRQ(ierr); 1107 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1108 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1109 tmatd = (Mat_SeqDense*)tmat->data; 1110 v = mat->v; v2 = tmatd->v; 1111 for (j=0; j<n; j++) { 1112 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1113 } 1114 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1115 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1116 *matout = tmat; 1117 } 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "MatEqual_SeqDense" 1123 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1124 { 1125 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1126 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1127 PetscInt i,j; 1128 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1129 1130 PetscFunctionBegin; 1131 if (A1->rmap.n != A2->rmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1132 if (A1->cmap.n != A2->cmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1133 for (i=0; i<A1->rmap.n; i++) { 1134 v1 = mat1->v+i; v2 = mat2->v+i; 1135 for (j=0; j<A1->cmap.n; j++) { 1136 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1137 v1 += mat1->lda; v2 += mat2->lda; 1138 } 1139 } 1140 *flg = PETSC_TRUE; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1146 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1147 { 1148 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1149 PetscErrorCode ierr; 1150 PetscInt i,n,len; 1151 PetscScalar *x,zero = 0.0; 1152 1153 PetscFunctionBegin; 1154 ierr = VecSet(v,zero);CHKERRQ(ierr); 1155 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1156 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1157 len = PetscMin(A->rmap.n,A->cmap.n); 1158 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1159 for (i=0; i<len; i++) { 1160 x[i] = mat->v[i*mat->lda + i]; 1161 } 1162 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1168 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1169 { 1170 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1171 PetscScalar *l,*r,x,*v; 1172 PetscErrorCode ierr; 1173 PetscInt i,j,m = A->rmap.n,n = A->cmap.n; 1174 1175 PetscFunctionBegin; 1176 if (ll) { 1177 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1178 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1179 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1180 for (i=0; i<m; i++) { 1181 x = l[i]; 1182 v = mat->v + i; 1183 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1184 } 1185 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1186 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1187 } 1188 if (rr) { 1189 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1190 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1191 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1192 for (i=0; i<n; i++) { 1193 x = r[i]; 1194 v = mat->v + i*m; 1195 for (j=0; j<m; j++) { (*v++) *= x;} 1196 } 1197 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1198 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1199 } 1200 PetscFunctionReturn(0); 1201 } 1202 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "MatNorm_SeqDense" 1205 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1206 { 1207 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1208 PetscScalar *v = mat->v; 1209 PetscReal sum = 0.0; 1210 PetscInt lda=mat->lda,m=A->rmap.n,i,j; 1211 PetscErrorCode ierr; 1212 1213 PetscFunctionBegin; 1214 if (type == NORM_FROBENIUS) { 1215 if (lda>m) { 1216 for (j=0; j<A->cmap.n; j++) { 1217 v = mat->v+j*lda; 1218 for (i=0; i<m; i++) { 1219 #if defined(PETSC_USE_COMPLEX) 1220 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1221 #else 1222 sum += (*v)*(*v); v++; 1223 #endif 1224 } 1225 } 1226 } else { 1227 for (i=0; i<A->cmap.n*A->rmap.n; i++) { 1228 #if defined(PETSC_USE_COMPLEX) 1229 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1230 #else 1231 sum += (*v)*(*v); v++; 1232 #endif 1233 } 1234 } 1235 *nrm = sqrt(sum); 1236 ierr = PetscLogFlops(2*A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1237 } else if (type == NORM_1) { 1238 *nrm = 0.0; 1239 for (j=0; j<A->cmap.n; j++) { 1240 v = mat->v + j*mat->lda; 1241 sum = 0.0; 1242 for (i=0; i<A->rmap.n; i++) { 1243 sum += PetscAbsScalar(*v); v++; 1244 } 1245 if (sum > *nrm) *nrm = sum; 1246 } 1247 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1248 } else if (type == NORM_INFINITY) { 1249 *nrm = 0.0; 1250 for (j=0; j<A->rmap.n; j++) { 1251 v = mat->v + j; 1252 sum = 0.0; 1253 for (i=0; i<A->cmap.n; i++) { 1254 sum += PetscAbsScalar(*v); v += mat->lda; 1255 } 1256 if (sum > *nrm) *nrm = sum; 1257 } 1258 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1259 } else { 1260 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1261 } 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatSetOption_SeqDense" 1267 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1268 { 1269 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1270 PetscErrorCode ierr; 1271 1272 PetscFunctionBegin; 1273 switch (op) { 1274 case MAT_ROW_ORIENTED: 1275 aij->roworiented = flg; 1276 break; 1277 case MAT_NEW_NONZERO_LOCATIONS: 1278 case MAT_NEW_NONZERO_LOCATION_ERR: 1279 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1280 case MAT_NEW_DIAGONALS: 1281 case MAT_IGNORE_OFF_PROC_ENTRIES: 1282 case MAT_USE_HASH_TABLE: 1283 case MAT_SYMMETRIC: 1284 case MAT_STRUCTURALLY_SYMMETRIC: 1285 case MAT_HERMITIAN: 1286 case MAT_SYMMETRY_ETERNAL: 1287 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1288 break; 1289 default: 1290 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatZeroEntries_SeqDense" 1297 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1298 { 1299 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1300 PetscErrorCode ierr; 1301 PetscInt lda=l->lda,m=A->rmap.n,j; 1302 1303 PetscFunctionBegin; 1304 if (lda>m) { 1305 for (j=0; j<A->cmap.n; j++) { 1306 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1307 } 1308 } else { 1309 ierr = PetscMemzero(l->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1310 } 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "MatZeroRows_SeqDense" 1316 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1317 { 1318 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1319 PetscInt n = A->cmap.n,i,j; 1320 PetscScalar *slot; 1321 1322 PetscFunctionBegin; 1323 for (i=0; i<N; i++) { 1324 slot = l->v + rows[i]; 1325 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1326 } 1327 if (diag != 0.0) { 1328 for (i=0; i<N; i++) { 1329 slot = l->v + (n+1)*rows[i]; 1330 *slot = diag; 1331 } 1332 } 1333 PetscFunctionReturn(0); 1334 } 1335 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "MatGetArray_SeqDense" 1338 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1339 { 1340 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1341 1342 PetscFunctionBegin; 1343 if (mat->lda != A->rmap.n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1344 *array = mat->v; 1345 PetscFunctionReturn(0); 1346 } 1347 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "MatRestoreArray_SeqDense" 1350 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1351 { 1352 PetscFunctionBegin; 1353 *array = 0; /* user cannot accidently use the array later */ 1354 PetscFunctionReturn(0); 1355 } 1356 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1359 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1360 { 1361 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1362 PetscErrorCode ierr; 1363 PetscInt i,j,*irow,*icol,nrows,ncols; 1364 PetscScalar *av,*bv,*v = mat->v; 1365 Mat newmat; 1366 1367 PetscFunctionBegin; 1368 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1369 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1370 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1371 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1372 1373 /* Check submatrixcall */ 1374 if (scall == MAT_REUSE_MATRIX) { 1375 PetscInt n_cols,n_rows; 1376 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1377 if (n_rows != nrows || n_cols != ncols) { 1378 /* resize the result result matrix to match number of requested rows/columns */ 1379 ierr = MatSetSizes(*B,nrows,nrows,nrows,nrows);CHKERRQ(ierr); 1380 } 1381 newmat = *B; 1382 } else { 1383 /* Create and fill new matrix */ 1384 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1385 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1386 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1387 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1388 } 1389 1390 /* Now extract the data pointers and do the copy,column at a time */ 1391 bv = ((Mat_SeqDense*)newmat->data)->v; 1392 1393 for (i=0; i<ncols; i++) { 1394 av = v + mat->lda*icol[i]; 1395 for (j=0; j<nrows; j++) { 1396 *bv++ = av[irow[j]]; 1397 } 1398 } 1399 1400 /* Assemble the matrices so that the correct flags are set */ 1401 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1402 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1403 1404 /* Free work space */ 1405 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1406 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1407 *B = newmat; 1408 PetscFunctionReturn(0); 1409 } 1410 1411 #undef __FUNCT__ 1412 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1413 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1414 { 1415 PetscErrorCode ierr; 1416 PetscInt i; 1417 1418 PetscFunctionBegin; 1419 if (scall == MAT_INITIAL_MATRIX) { 1420 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1421 } 1422 1423 for (i=0; i<n; i++) { 1424 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1425 } 1426 PetscFunctionReturn(0); 1427 } 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1431 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1432 { 1433 PetscFunctionBegin; 1434 PetscFunctionReturn(0); 1435 } 1436 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1439 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1440 { 1441 PetscFunctionBegin; 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #undef __FUNCT__ 1446 #define __FUNCT__ "MatCopy_SeqDense" 1447 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1448 { 1449 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1450 PetscErrorCode ierr; 1451 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap.n,n=A->cmap.n, j; 1452 1453 PetscFunctionBegin; 1454 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1455 if (A->ops->copy != B->ops->copy) { 1456 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1457 PetscFunctionReturn(0); 1458 } 1459 if (m != B->rmap.n || n != B->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1460 if (lda1>m || lda2>m) { 1461 for (j=0; j<n; j++) { 1462 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1463 } 1464 } else { 1465 ierr = PetscMemcpy(b->v,a->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1466 } 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1472 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1473 { 1474 PetscErrorCode ierr; 1475 1476 PetscFunctionBegin; 1477 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "MatSetSizes_SeqDense" 1483 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1484 { 1485 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1486 PetscErrorCode ierr; 1487 PetscFunctionBegin; 1488 /* this will not be called before lda, Mmax, and Nmax have been set */ 1489 m = PetscMax(m,M); 1490 n = PetscMax(n,N); 1491 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); 1492 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); 1493 A->rmap.n = A->rmap.n = m; 1494 A->cmap.n = A->cmap.N = n; 1495 if (a->changelda) a->lda = m; 1496 ierr = PetscMemzero(a->v,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 1501 /* ----------------------------------------------------------------*/ 1502 #undef __FUNCT__ 1503 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1504 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1505 { 1506 PetscErrorCode ierr; 1507 1508 PetscFunctionBegin; 1509 if (scall == MAT_INITIAL_MATRIX){ 1510 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1511 } 1512 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1513 PetscFunctionReturn(0); 1514 } 1515 1516 #undef __FUNCT__ 1517 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1518 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1519 { 1520 PetscErrorCode ierr; 1521 PetscInt m=A->rmap.n,n=B->cmap.n; 1522 Mat Cmat; 1523 1524 PetscFunctionBegin; 1525 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); 1526 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1527 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1528 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1529 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1530 Cmat->assembled = PETSC_TRUE; 1531 *C = Cmat; 1532 PetscFunctionReturn(0); 1533 } 1534 1535 #undef __FUNCT__ 1536 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1537 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1538 { 1539 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1540 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1541 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1542 PetscBLASInt m,n,k; 1543 PetscScalar _DOne=1.0,_DZero=0.0; 1544 1545 PetscFunctionBegin; 1546 m = PetscBLASIntCast(A->rmap.n); 1547 n = PetscBLASIntCast(B->cmap.n); 1548 k = PetscBLASIntCast(A->cmap.n); 1549 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 #undef __FUNCT__ 1554 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1555 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1556 { 1557 PetscErrorCode ierr; 1558 1559 PetscFunctionBegin; 1560 if (scall == MAT_INITIAL_MATRIX){ 1561 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1562 } 1563 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1564 PetscFunctionReturn(0); 1565 } 1566 1567 #undef __FUNCT__ 1568 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1569 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1570 { 1571 PetscErrorCode ierr; 1572 PetscInt m=A->cmap.n,n=B->cmap.n; 1573 Mat Cmat; 1574 1575 PetscFunctionBegin; 1576 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); 1577 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1578 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1579 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1580 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1581 Cmat->assembled = PETSC_TRUE; 1582 *C = Cmat; 1583 PetscFunctionReturn(0); 1584 } 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1588 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1589 { 1590 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1591 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1592 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1593 PetscBLASInt m,n,k; 1594 PetscScalar _DOne=1.0,_DZero=0.0; 1595 1596 PetscFunctionBegin; 1597 m = PetscBLASIntCast(A->rmap.n); 1598 n = PetscBLASIntCast(B->cmap.n); 1599 k = PetscBLASIntCast(A->cmap.n); 1600 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1601 PetscFunctionReturn(0); 1602 } 1603 1604 #undef __FUNCT__ 1605 #define __FUNCT__ "MatGetRowMax_SeqDense" 1606 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1607 { 1608 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1609 PetscErrorCode ierr; 1610 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1611 PetscScalar *x; 1612 MatScalar *aa = a->v; 1613 1614 PetscFunctionBegin; 1615 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1616 1617 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1618 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1619 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1620 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1621 for (i=0; i<m; i++) { 1622 x[i] = aa[i]; if (idx) idx[i] = 0; 1623 for (j=1; j<n; j++){ 1624 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1625 } 1626 } 1627 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1633 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1634 { 1635 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1636 PetscErrorCode ierr; 1637 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1638 PetscScalar *x; 1639 PetscReal atmp; 1640 MatScalar *aa = a->v; 1641 1642 PetscFunctionBegin; 1643 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1644 1645 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1646 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1647 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1648 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1649 for (i=0; i<m; i++) { 1650 x[i] = PetscAbsScalar(aa[i]); if (idx) idx[i] = 0; 1651 for (j=1; j<n; j++){ 1652 atmp = PetscAbsScalar(aa[i+m*j]); 1653 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1654 } 1655 } 1656 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 #undef __FUNCT__ 1661 #define __FUNCT__ "MatGetRowMin_SeqDense" 1662 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1663 { 1664 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1665 PetscErrorCode ierr; 1666 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1667 PetscScalar *x; 1668 MatScalar *aa = a->v; 1669 1670 PetscFunctionBegin; 1671 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1672 1673 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1674 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1675 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1676 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1677 for (i=0; i<m; i++) { 1678 x[i] = aa[i]; if (idx) idx[i] = 0; 1679 for (j=1; j<n; j++){ 1680 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1681 } 1682 } 1683 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1684 PetscFunctionReturn(0); 1685 } 1686 1687 #undef __FUNCT__ 1688 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1689 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1690 { 1691 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1692 PetscErrorCode ierr; 1693 PetscScalar *x; 1694 1695 PetscFunctionBegin; 1696 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1697 1698 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1699 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1700 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1701 PetscFunctionReturn(0); 1702 } 1703 1704 /* -------------------------------------------------------------------*/ 1705 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1706 MatGetRow_SeqDense, 1707 MatRestoreRow_SeqDense, 1708 MatMult_SeqDense, 1709 /* 4*/ MatMultAdd_SeqDense, 1710 MatMultTranspose_SeqDense, 1711 MatMultTransposeAdd_SeqDense, 1712 MatSolve_SeqDense, 1713 MatSolveAdd_SeqDense, 1714 MatSolveTranspose_SeqDense, 1715 /*10*/ MatSolveTransposeAdd_SeqDense, 1716 MatLUFactor_SeqDense, 1717 MatCholeskyFactor_SeqDense, 1718 MatRelax_SeqDense, 1719 MatTranspose_SeqDense, 1720 /*15*/ MatGetInfo_SeqDense, 1721 MatEqual_SeqDense, 1722 MatGetDiagonal_SeqDense, 1723 MatDiagonalScale_SeqDense, 1724 MatNorm_SeqDense, 1725 /*20*/ MatAssemblyBegin_SeqDense, 1726 MatAssemblyEnd_SeqDense, 1727 0, 1728 MatSetOption_SeqDense, 1729 MatZeroEntries_SeqDense, 1730 /*25*/ MatZeroRows_SeqDense, 1731 MatLUFactorSymbolic_SeqDense, 1732 MatLUFactorNumeric_SeqDense, 1733 MatCholeskyFactorSymbolic_SeqDense, 1734 MatCholeskyFactorNumeric_SeqDense, 1735 /*30*/ MatSetUpPreallocation_SeqDense, 1736 0, 1737 0, 1738 MatGetArray_SeqDense, 1739 MatRestoreArray_SeqDense, 1740 /*35*/ MatDuplicate_SeqDense, 1741 0, 1742 0, 1743 0, 1744 0, 1745 /*40*/ MatAXPY_SeqDense, 1746 MatGetSubMatrices_SeqDense, 1747 0, 1748 MatGetValues_SeqDense, 1749 MatCopy_SeqDense, 1750 /*45*/ MatGetRowMax_SeqDense, 1751 MatScale_SeqDense, 1752 0, 1753 0, 1754 0, 1755 /*50*/ 0, 1756 0, 1757 0, 1758 0, 1759 0, 1760 /*55*/ 0, 1761 0, 1762 0, 1763 0, 1764 0, 1765 /*60*/ 0, 1766 MatDestroy_SeqDense, 1767 MatView_SeqDense, 1768 0, 1769 0, 1770 /*65*/ 0, 1771 0, 1772 0, 1773 0, 1774 0, 1775 /*70*/ MatGetRowMaxAbs_SeqDense, 1776 0, 1777 0, 1778 0, 1779 0, 1780 /*75*/ 0, 1781 0, 1782 0, 1783 0, 1784 0, 1785 /*80*/ 0, 1786 0, 1787 0, 1788 0, 1789 /*84*/ MatLoad_SeqDense, 1790 0, 1791 MatIsHermitian_SeqDense, 1792 0, 1793 0, 1794 0, 1795 /*90*/ MatMatMult_SeqDense_SeqDense, 1796 MatMatMultSymbolic_SeqDense_SeqDense, 1797 MatMatMultNumeric_SeqDense_SeqDense, 1798 0, 1799 0, 1800 /*95*/ 0, 1801 MatMatMultTranspose_SeqDense_SeqDense, 1802 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1803 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1804 0, 1805 /*100*/0, 1806 0, 1807 0, 1808 0, 1809 MatSetSizes_SeqDense, 1810 0, 1811 0, 1812 0, 1813 0, 1814 0, 1815 /*110*/0, 1816 0, 1817 MatGetRowMin_SeqDense, 1818 MatGetColumnVector_SeqDense 1819 }; 1820 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "MatCreateSeqDense" 1823 /*@C 1824 MatCreateSeqDense - Creates a sequential dense matrix that 1825 is stored in column major order (the usual Fortran 77 manner). Many 1826 of the matrix operations use the BLAS and LAPACK routines. 1827 1828 Collective on MPI_Comm 1829 1830 Input Parameters: 1831 + comm - MPI communicator, set to PETSC_COMM_SELF 1832 . m - number of rows 1833 . n - number of columns 1834 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1835 to control all matrix memory allocation. 1836 1837 Output Parameter: 1838 . A - the matrix 1839 1840 Notes: 1841 The data input variable is intended primarily for Fortran programmers 1842 who wish to allocate their own matrix memory space. Most users should 1843 set data=PETSC_NULL. 1844 1845 Level: intermediate 1846 1847 .keywords: dense, matrix, LAPACK, BLAS 1848 1849 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1850 @*/ 1851 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1852 { 1853 PetscErrorCode ierr; 1854 1855 PetscFunctionBegin; 1856 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1857 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1858 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1859 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1860 PetscFunctionReturn(0); 1861 } 1862 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1865 /*@C 1866 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1867 1868 Collective on MPI_Comm 1869 1870 Input Parameters: 1871 + A - the matrix 1872 - data - the array (or PETSC_NULL) 1873 1874 Notes: 1875 The data input variable is intended primarily for Fortran programmers 1876 who wish to allocate their own matrix memory space. Most users should 1877 need not call this routine. 1878 1879 Level: intermediate 1880 1881 .keywords: dense, matrix, LAPACK, BLAS 1882 1883 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1884 @*/ 1885 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1886 { 1887 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1888 1889 PetscFunctionBegin; 1890 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1891 if (f) { 1892 ierr = (*f)(B,data);CHKERRQ(ierr); 1893 } 1894 PetscFunctionReturn(0); 1895 } 1896 1897 EXTERN_C_BEGIN 1898 #undef __FUNCT__ 1899 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1900 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1901 { 1902 Mat_SeqDense *b; 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 B->preallocated = PETSC_TRUE; 1907 b = (Mat_SeqDense*)B->data; 1908 if (!data) { /* petsc-allocated storage */ 1909 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1910 ierr = PetscMalloc((b->lda*b->Nmax+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1911 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1912 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1913 b->user_alloc = PETSC_FALSE; 1914 } else { /* user-allocated storage */ 1915 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1916 b->v = data; 1917 b->user_alloc = PETSC_TRUE; 1918 } 1919 PetscFunctionReturn(0); 1920 } 1921 EXTERN_C_END 1922 1923 #undef __FUNCT__ 1924 #define __FUNCT__ "MatSeqDenseSetLDA" 1925 /*@C 1926 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1927 1928 Input parameter: 1929 + A - the matrix 1930 - lda - the leading dimension 1931 1932 Notes: 1933 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1934 it asserts that the preallocation has a leading dimension (the LDA parameter 1935 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1936 1937 Level: intermediate 1938 1939 .keywords: dense, matrix, LAPACK, BLAS 1940 1941 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 1942 @*/ 1943 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 1944 { 1945 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1946 1947 PetscFunctionBegin; 1948 if (lda < B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap.n); 1949 b->lda = lda; 1950 b->changelda = PETSC_FALSE; 1951 b->Mmax = PetscMax(b->Mmax,lda); 1952 PetscFunctionReturn(0); 1953 } 1954 1955 /*MC 1956 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1957 1958 Options Database Keys: 1959 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1960 1961 Level: beginner 1962 1963 .seealso: MatCreateSeqDense() 1964 1965 M*/ 1966 1967 EXTERN_C_BEGIN 1968 #undef __FUNCT__ 1969 #define __FUNCT__ "MatCreate_SeqDense" 1970 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 1971 { 1972 Mat_SeqDense *b; 1973 PetscErrorCode ierr; 1974 PetscMPIInt size; 1975 1976 PetscFunctionBegin; 1977 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1978 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1979 1980 B->rmap.bs = B->cmap.bs = 1; 1981 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 1982 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 1983 1984 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 1985 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1986 B->factor = 0; 1987 B->mapping = 0; 1988 B->data = (void*)b; 1989 1990 1991 b->pivots = 0; 1992 b->roworiented = PETSC_TRUE; 1993 b->v = 0; 1994 b->lda = B->rmap.n; 1995 b->changelda = PETSC_FALSE; 1996 b->Mmax = B->rmap.n; 1997 b->Nmax = B->cmap.n; 1998 1999 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2000 "MatSeqDenseSetPreallocation_SeqDense", 2001 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2002 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2003 "MatMatMult_SeqAIJ_SeqDense", 2004 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2005 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2006 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2007 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2009 "MatMatMultNumeric_SeqAIJ_SeqDense", 2010 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2011 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2012 PetscFunctionReturn(0); 2013 } 2014 2015 2016 EXTERN_C_END 2017