1 /*$Id: dense.c,v 1.208 2001/09/07 20:09:16 bsmith Exp $*/ 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include "src/mat/impls/dense/seq/dense.h" 7 #include "petscblaslapack.h" 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatAXPY_SeqDense" 11 int MatAXPY_SeqDense(const PetscScalar *alpha,Mat X,Mat Y,MatStructure str) 12 { 13 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 14 int N = X->m*X->n,m=X->m,ldax=x->lda,lday=y->lda, j,one = 1; 15 16 PetscFunctionBegin; 17 if (X->m != Y->m || X->n != Y->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 18 if (ldax>m || lday>m) { 19 for (j=0; j<X->n; j++) { 20 BLaxpy_(&m,(PetscScalar*)alpha,x->v+j*ldax,&one,y->v+j*lday,&one); 21 } 22 } else { 23 BLaxpy_(&N,(PetscScalar*)alpha,x->v,&one,y->v,&one); 24 } 25 PetscLogFlops(2*N-1); 26 PetscFunctionReturn(0); 27 } 28 29 #undef __FUNCT__ 30 #define __FUNCT__ "MatGetInfo_SeqDense" 31 int MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 32 { 33 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 34 int i,N = A->m*A->n,count = 0; 35 PetscScalar *v = a->v; 36 37 PetscFunctionBegin; 38 for (i=0; i<N; i++) {if (*v != 0.0) count++; v++;} 39 40 info->rows_global = (double)A->m; 41 info->columns_global = (double)A->n; 42 info->rows_local = (double)A->m; 43 info->columns_local = (double)A->n; 44 info->block_size = 1.0; 45 info->nz_allocated = (double)N; 46 info->nz_used = (double)count; 47 info->nz_unneeded = (double)(N-count); 48 info->assemblies = (double)A->num_ass; 49 info->mallocs = 0; 50 info->memory = A->mem; 51 info->fill_ratio_given = 0; 52 info->fill_ratio_needed = 0; 53 info->factor_mallocs = 0; 54 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "MatScale_SeqDense" 60 int MatScale_SeqDense(const PetscScalar *alpha,Mat A) 61 { 62 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 63 int one = 1,lda = a->lda,j,nz; 64 65 PetscFunctionBegin; 66 if (lda>A->m) { 67 nz = A->m; 68 for (j=0; j<A->n; j++) { 69 BLscal_(&nz,(PetscScalar*)alpha,a->v+j*lda,&one); 70 } 71 } else { 72 nz = A->m*A->n; 73 BLscal_(&nz,(PetscScalar*)alpha,a->v,&one); 74 } 75 PetscLogFlops(nz); 76 PetscFunctionReturn(0); 77 } 78 79 /* ---------------------------------------------------------------*/ 80 /* COMMENT: I have chosen to hide row permutation in the pivots, 81 rather than put it in the Mat->row slot.*/ 82 #undef __FUNCT__ 83 #define __FUNCT__ "MatLUFactor_SeqDense" 84 int MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo) 85 { 86 #if defined(PETSC_MISSING_LAPACK_GETRF) 87 PetscFunctionBegin; 88 SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 89 #else 90 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 91 int info,ierr; 92 93 PetscFunctionBegin; 94 if (!mat->pivots) { 95 ierr = PetscMalloc((A->m+1)*sizeof(int),&mat->pivots);CHKERRQ(ierr); 96 PetscLogObjectMemory(A,A->m*sizeof(int)); 97 } 98 A->factor = FACTOR_LU; 99 if (!A->m || !A->n) PetscFunctionReturn(0); 100 LAgetrf_(&A->m,&A->n,mat->v,&mat->lda,mat->pivots,&info); 101 if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization"); 102 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 103 PetscLogFlops((2*A->n*A->n*A->n)/3); 104 #endif 105 PetscFunctionReturn(0); 106 } 107 108 #undef __FUNCT__ 109 #define __FUNCT__ "MatDuplicate_SeqDense" 110 int MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 111 { 112 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 113 int lda = mat->lda,j,m,ierr; 114 Mat newi; 115 116 PetscFunctionBegin; 117 ierr = MatCreateSeqDense(A->comm,A->m,A->n,PETSC_NULL,&newi);CHKERRQ(ierr); 118 if (cpvalues == MAT_COPY_VALUES) { 119 l = (Mat_SeqDense*)newi->data; 120 if (lda>A->m) { 121 m = A->m; 122 for (j=0; j<A->n; j++) { 123 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 124 } 125 } else { 126 ierr = PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 127 } 128 } 129 newi->assembled = PETSC_TRUE; 130 *newmat = newi; 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 136 int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatFactorInfo *info,Mat *fact) 137 { 138 int ierr; 139 140 PetscFunctionBegin; 141 ierr = MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 147 int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 148 { 149 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l = (Mat_SeqDense*)(*fact)->data; 150 int lda1=mat->lda,lda2=l->lda, m=A->m,n=A->n, j,ierr; 151 152 PetscFunctionBegin; 153 /* copy the numerical values */ 154 if (lda1>m || lda2>m ) { 155 for (j=0; j<n; j++) { 156 ierr = PetscMemcpy(l->v+j*lda2,mat->v+j*lda1,m*sizeof(PetscScalar)); CHKERRQ(ierr); 157 } 158 } else { 159 ierr = PetscMemcpy(l->v,mat->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 160 } 161 (*fact)->factor = 0; 162 ierr = MatLUFactor(*fact,0,0,PETSC_NULL);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 168 int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,MatFactorInfo *info,Mat *fact) 169 { 170 int ierr; 171 172 PetscFunctionBegin; 173 ierr = MatConvert(A,MATSAME,fact);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 179 int MatCholeskyFactor_SeqDense(Mat A,IS perm,MatFactorInfo *factinfo) 180 { 181 #if defined(PETSC_MISSING_LAPACK_POTRF) 182 PetscFunctionBegin; 183 SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 184 #else 185 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 186 int info,ierr; 187 188 PetscFunctionBegin; 189 if (mat->pivots) { 190 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 191 PetscLogObjectMemory(A,-A->m*sizeof(int)); 192 mat->pivots = 0; 193 } 194 if (!A->m || !A->n) PetscFunctionReturn(0); 195 LApotrf_("L",&A->n,mat->v,&mat->lda,&info); 196 if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %d",info-1); 197 A->factor = FACTOR_CHOLESKY; 198 PetscLogFlops((A->n*A->n*A->n)/3); 199 #endif 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 205 int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 206 { 207 int ierr; 208 MatFactorInfo info; 209 210 PetscFunctionBegin; 211 info.fill = 1.0; 212 ierr = MatCholeskyFactor_SeqDense(*fact,0,&info);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "MatSolve_SeqDense" 218 int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 219 { 220 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 221 int one = 1,info,ierr; 222 PetscScalar *x,*y; 223 224 PetscFunctionBegin; 225 if (!A->m || !A->n) PetscFunctionReturn(0); 226 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 227 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 228 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 229 if (A->factor == FACTOR_LU) { 230 #if defined(PETSC_MISSING_LAPACK_GETRS) 231 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 232 #else 233 LAgetrs_("N",&A->m,&one,mat->v,&mat->lda,mat->pivots,y,&A->m,&info); 234 if (info) SETERRQ(PETSC_ERR_LIB,"MBad solve"); 235 #endif 236 } else if (A->factor == FACTOR_CHOLESKY){ 237 #if defined(PETSC_MISSING_LAPACK_POTRS) 238 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 239 #else 240 LApotrs_("L",&A->m,&one,mat->v,&mat->lda,y,&A->m,&info); 241 if (info) SETERRQ(PETSC_ERR_LIB,"MBad solve"); 242 #endif 243 } 244 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 245 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 246 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 247 PetscLogFlops(2*A->n*A->n - A->n); 248 PetscFunctionReturn(0); 249 } 250 251 #undef __FUNCT__ 252 #define __FUNCT__ "MatSolveTranspose_SeqDense" 253 int MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 254 { 255 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 256 int ierr,one = 1,info; 257 PetscScalar *x,*y; 258 259 PetscFunctionBegin; 260 if (!A->m || !A->n) PetscFunctionReturn(0); 261 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 262 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 263 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 264 /* assume if pivots exist then use LU; else Cholesky */ 265 if (mat->pivots) { 266 #if defined(PETSC_MISSING_LAPACK_GETRS) 267 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 268 #else 269 LAgetrs_("T",&A->m,&one,mat->v,&mat->lda,mat->pivots,y,&A->m,&info); 270 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 271 #endif 272 } else { 273 #if defined(PETSC_MISSING_LAPACK_POTRS) 274 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 275 #else 276 LApotrs_("L",&A->m,&one,mat->v,&mat->lda,y,&A->m,&info); 277 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 278 #endif 279 } 280 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 281 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 282 PetscLogFlops(2*A->n*A->n - A->n); 283 PetscFunctionReturn(0); 284 } 285 286 #undef __FUNCT__ 287 #define __FUNCT__ "MatSolveAdd_SeqDense" 288 int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 289 { 290 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 291 int ierr,one = 1,info; 292 PetscScalar *x,*y,sone = 1.0; 293 Vec tmp = 0; 294 295 PetscFunctionBegin; 296 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 297 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 298 if (!A->m || !A->n) PetscFunctionReturn(0); 299 if (yy == zz) { 300 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 301 PetscLogObjectParent(A,tmp); 302 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 303 } 304 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 305 /* assume if pivots exist then use LU; else Cholesky */ 306 if (mat->pivots) { 307 #if defined(PETSC_MISSING_LAPACK_GETRS) 308 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 309 #else 310 LAgetrs_("N",&A->m,&one,mat->v,&mat->lda,mat->pivots,y,&A->m,&info); 311 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 312 #endif 313 } else { 314 #if defined(PETSC_MISSING_LAPACK_POTRS) 315 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 316 #else 317 LApotrs_("L",&A->m,&one,mat->v,&mat->lda,y,&A->m,&info); 318 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 319 #endif 320 } 321 if (tmp) {ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 322 else {ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr);} 323 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 324 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 325 PetscLogFlops(2*A->n*A->n); 326 PetscFunctionReturn(0); 327 } 328 329 #undef __FUNCT__ 330 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 331 int MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 332 { 333 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 334 int one = 1,info,ierr; 335 PetscScalar *x,*y,sone = 1.0; 336 Vec tmp; 337 338 PetscFunctionBegin; 339 if (!A->m || !A->n) PetscFunctionReturn(0); 340 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 341 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 342 if (yy == zz) { 343 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 344 PetscLogObjectParent(A,tmp); 345 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 346 } 347 ierr = PetscMemcpy(y,x,A->m*sizeof(PetscScalar));CHKERRQ(ierr); 348 /* assume if pivots exist then use LU; else Cholesky */ 349 if (mat->pivots) { 350 #if defined(PETSC_MISSING_LAPACK_GETRS) 351 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 352 #else 353 LAgetrs_("T",&A->m,&one,mat->v,&mat->lda,mat->pivots,y,&A->m,&info); 354 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 355 #endif 356 } else { 357 #if defined(PETSC_MISSING_LAPACK_POTRS) 358 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 359 #else 360 LApotrs_("L",&A->m,&one,mat->v,&mat->lda,y,&A->m,&info); 361 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 362 #endif 363 } 364 if (tmp) { 365 ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); 366 ierr = VecDestroy(tmp);CHKERRQ(ierr); 367 } else { 368 ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr); 369 } 370 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 371 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 372 PetscLogFlops(2*A->n*A->n); 373 PetscFunctionReturn(0); 374 } 375 /* ------------------------------------------------------------------*/ 376 #undef __FUNCT__ 377 #define __FUNCT__ "MatRelax_SeqDense" 378 int MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag, 379 PetscReal shift,int its,int lits,Vec xx) 380 { 381 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 382 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 383 int ierr,m = A->m,i; 384 #if !defined(PETSC_USE_COMPLEX) 385 int o = 1; 386 #endif 387 388 PetscFunctionBegin; 389 if (flag & SOR_ZERO_INITIAL_GUESS) { 390 /* this is a hack fix, should have another version without the second BLdot */ 391 ierr = VecSet(&zero,xx);CHKERRQ(ierr); 392 } 393 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 394 ierr = VecGetArrayFast(bb,&b);CHKERRQ(ierr); 395 its = its*lits; 396 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 397 while (its--) { 398 if (flag & SOR_FORWARD_SWEEP){ 399 for (i=0; i<m; i++) { 400 #if defined(PETSC_USE_COMPLEX) 401 /* cannot use BLAS dot for complex because compiler/linker is 402 not happy about returning a double complex */ 403 int _i; 404 PetscScalar sum = b[i]; 405 for (_i=0; _i<m; _i++) { 406 sum -= PetscConj(v[i+_i*m])*x[_i]; 407 } 408 xt = sum; 409 #else 410 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 411 #endif 412 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 413 } 414 } 415 if (flag & SOR_BACKWARD_SWEEP) { 416 for (i=m-1; i>=0; i--) { 417 #if defined(PETSC_USE_COMPLEX) 418 /* cannot use BLAS dot for complex because compiler/linker is 419 not happy about returning a double complex */ 420 int _i; 421 PetscScalar sum = b[i]; 422 for (_i=0; _i<m; _i++) { 423 sum -= PetscConj(v[i+_i*m])*x[_i]; 424 } 425 xt = sum; 426 #else 427 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 428 #endif 429 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 430 } 431 } 432 } 433 ierr = VecRestoreArrayFast(bb,&b);CHKERRQ(ierr); 434 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 435 PetscFunctionReturn(0); 436 } 437 438 /* -----------------------------------------------------------------*/ 439 #undef __FUNCT__ 440 #define __FUNCT__ "MatMultTranspose_SeqDense" 441 int MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 442 { 443 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 444 PetscScalar *v = mat->v,*x,*y; 445 int ierr,_One=1; 446 PetscScalar _DOne=1.0,_DZero=0.0; 447 448 PetscFunctionBegin; 449 if (!A->m || !A->n) PetscFunctionReturn(0); 450 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 451 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 452 LAgemv_("T",&(A->m),&(A->n),&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 453 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 454 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 455 PetscLogFlops(2*A->m*A->n - A->n); 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "MatMult_SeqDense" 461 int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 462 { 463 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 464 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 465 int ierr,_One=1; 466 467 PetscFunctionBegin; 468 if (!A->m || !A->n) PetscFunctionReturn(0); 469 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 470 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 471 LAgemv_("N",&(A->m),&(A->n),&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 472 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 473 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 474 PetscLogFlops(2*A->m*A->n - A->m); 475 PetscFunctionReturn(0); 476 } 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "MatMultAdd_SeqDense" 480 int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 481 { 482 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 483 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 484 int ierr,_One=1; 485 486 PetscFunctionBegin; 487 if (!A->m || !A->n) PetscFunctionReturn(0); 488 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 489 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 490 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 491 LAgemv_("N",&(A->m),&(A->n),&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 492 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 493 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 494 PetscLogFlops(2*A->m*A->n); 495 PetscFunctionReturn(0); 496 } 497 498 #undef __FUNCT__ 499 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 500 int MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 501 { 502 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 503 PetscScalar *v = mat->v,*x,*y; 504 int ierr,_One=1; 505 PetscScalar _DOne=1.0; 506 507 PetscFunctionBegin; 508 if (!A->m || !A->n) PetscFunctionReturn(0); 509 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 510 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 511 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 512 LAgemv_("T",&(A->m),&(A->n),&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 513 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 514 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 515 PetscLogFlops(2*A->m*A->n); 516 PetscFunctionReturn(0); 517 } 518 519 /* -----------------------------------------------------------------*/ 520 #undef __FUNCT__ 521 #define __FUNCT__ "MatGetRow_SeqDense" 522 int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals) 523 { 524 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 525 PetscScalar *v; 526 int i,ierr; 527 528 PetscFunctionBegin; 529 *ncols = A->n; 530 if (cols) { 531 ierr = PetscMalloc((A->n+1)*sizeof(int),cols);CHKERRQ(ierr); 532 for (i=0; i<A->n; i++) (*cols)[i] = i; 533 } 534 if (vals) { 535 ierr = PetscMalloc((A->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 536 v = mat->v + row; 537 for (i=0; i<A->n; i++) {(*vals)[i] = *v; v += mat->lda;} 538 } 539 PetscFunctionReturn(0); 540 } 541 542 #undef __FUNCT__ 543 #define __FUNCT__ "MatRestoreRow_SeqDense" 544 int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,PetscScalar **vals) 545 { 546 int ierr; 547 PetscFunctionBegin; 548 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 549 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 550 PetscFunctionReturn(0); 551 } 552 /* ----------------------------------------------------------------*/ 553 #undef __FUNCT__ 554 #define __FUNCT__ "MatSetValues_SeqDense" 555 int MatSetValues_SeqDense(Mat A,int m,const int indexm[],int n,const int indexn[],const PetscScalar v[],InsertMode addv) 556 { 557 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 558 int i,j,idx=0; 559 560 PetscFunctionBegin; 561 if (!mat->roworiented) { 562 if (addv == INSERT_VALUES) { 563 for (j=0; j<n; j++) { 564 if (indexn[j] < 0) {idx += m; continue;} 565 #if defined(PETSC_USE_BOPT_g) 566 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",indexn[j],A->n-1); 567 #endif 568 for (i=0; i<m; i++) { 569 if (indexm[i] < 0) {idx++; continue;} 570 #if defined(PETSC_USE_BOPT_g) 571 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",indexm[i],A->m-1); 572 #endif 573 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 574 } 575 } 576 } else { 577 for (j=0; j<n; j++) { 578 if (indexn[j] < 0) {idx += m; continue;} 579 #if defined(PETSC_USE_BOPT_g) 580 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",indexn[j],A->n-1); 581 #endif 582 for (i=0; i<m; i++) { 583 if (indexm[i] < 0) {idx++; continue;} 584 #if defined(PETSC_USE_BOPT_g) 585 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",indexm[i],A->m-1); 586 #endif 587 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 588 } 589 } 590 } 591 } else { 592 if (addv == INSERT_VALUES) { 593 for (i=0; i<m; i++) { 594 if (indexm[i] < 0) { idx += n; continue;} 595 #if defined(PETSC_USE_BOPT_g) 596 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",indexm[j],A->m-1); 597 #endif 598 for (j=0; j<n; j++) { 599 if (indexn[j] < 0) { idx++; continue;} 600 #if defined(PETSC_USE_BOPT_g) 601 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",indexn[i],A->n-1); 602 #endif 603 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 604 } 605 } 606 } else { 607 for (i=0; i<m; i++) { 608 if (indexm[i] < 0) { idx += n; continue;} 609 #if defined(PETSC_USE_BOPT_g) 610 if (indexm[i] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",indexm[j],A->m-1); 611 #endif 612 for (j=0; j<n; j++) { 613 if (indexn[j] < 0) { idx++; continue;} 614 #if defined(PETSC_USE_BOPT_g) 615 if (indexn[j] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",indexn[i],A->n-1); 616 #endif 617 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 618 } 619 } 620 } 621 } 622 PetscFunctionReturn(0); 623 } 624 625 #undef __FUNCT__ 626 #define __FUNCT__ "MatGetValues_SeqDense" 627 int MatGetValues_SeqDense(Mat A,int m,const int indexm[],int n,const int indexn[],PetscScalar v[]) 628 { 629 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 630 int i,j; 631 PetscScalar *vpt = v; 632 633 PetscFunctionBegin; 634 /* row-oriented output */ 635 for (i=0; i<m; i++) { 636 for (j=0; j<n; j++) { 637 *vpt++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 638 } 639 } 640 PetscFunctionReturn(0); 641 } 642 643 /* -----------------------------------------------------------------*/ 644 645 #include "petscsys.h" 646 647 #undef __FUNCT__ 648 #define __FUNCT__ "MatLoad_SeqDense" 649 int MatLoad_SeqDense(PetscViewer viewer,const MatType type,Mat *A) 650 { 651 Mat_SeqDense *a; 652 Mat B; 653 int *scols,i,j,nz,ierr,fd,header[4],size; 654 int *rowlengths = 0,M,N,*cols; 655 PetscScalar *vals,*svals,*v,*w; 656 MPI_Comm comm = ((PetscObject)viewer)->comm; 657 658 PetscFunctionBegin; 659 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 660 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 661 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 662 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 663 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 664 M = header[1]; N = header[2]; nz = header[3]; 665 666 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 667 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 668 B = *A; 669 a = (Mat_SeqDense*)B->data; 670 v = a->v; 671 /* Allocate some temp space to read in the values and then flip them 672 from row major to column major */ 673 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 674 /* read in nonzero values */ 675 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 676 /* now flip the values and store them in the matrix*/ 677 for (j=0; j<N; j++) { 678 for (i=0; i<M; i++) { 679 *v++ =w[i*N+j]; 680 } 681 } 682 ierr = PetscFree(w);CHKERRQ(ierr); 683 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 684 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 685 } else { 686 /* read row lengths */ 687 ierr = PetscMalloc((M+1)*sizeof(int),&rowlengths);CHKERRQ(ierr); 688 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 689 690 /* create our matrix */ 691 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 692 B = *A; 693 a = (Mat_SeqDense*)B->data; 694 v = a->v; 695 696 /* read column indices and nonzeros */ 697 ierr = PetscMalloc((nz+1)*sizeof(int),&scols);CHKERRQ(ierr); 698 cols = scols; 699 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 700 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 701 vals = svals; 702 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 703 704 /* insert into matrix */ 705 for (i=0; i<M; i++) { 706 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 707 svals += rowlengths[i]; scols += rowlengths[i]; 708 } 709 ierr = PetscFree(vals);CHKERRQ(ierr); 710 ierr = PetscFree(cols);CHKERRQ(ierr); 711 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 712 713 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 714 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 715 } 716 PetscFunctionReturn(0); 717 } 718 719 #include "petscsys.h" 720 721 #undef __FUNCT__ 722 #define __FUNCT__ "MatView_SeqDense_ASCII" 723 static int MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 724 { 725 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 726 int ierr,i,j; 727 char *name; 728 PetscScalar *v; 729 PetscViewerFormat format; 730 731 PetscFunctionBegin; 732 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 733 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 734 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 735 PetscFunctionReturn(0); /* do nothing for now */ 736 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 737 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 738 for (i=0; i<A->m; i++) { 739 v = a->v + i; 740 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 741 for (j=0; j<A->n; j++) { 742 #if defined(PETSC_USE_COMPLEX) 743 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 744 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 745 } else if (PetscRealPart(*v)) { 746 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",j,PetscRealPart(*v));CHKERRQ(ierr); 747 } 748 #else 749 if (*v) { 750 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",j,*v);CHKERRQ(ierr); 751 } 752 #endif 753 v += a->lda; 754 } 755 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 756 } 757 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 758 } else { 759 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 760 #if defined(PETSC_USE_COMPLEX) 761 PetscTruth allreal = PETSC_TRUE; 762 /* determine if matrix has all real values */ 763 v = a->v; 764 for (i=0; i<A->m*A->n; i++) { 765 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 766 } 767 #endif 768 if (format == PETSC_VIEWER_ASCII_MATLAB) { 769 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 770 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",A->m,A->n);CHKERRQ(ierr); 771 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%d,%d);\n",name,A->m,A->n);CHKERRQ(ierr); 772 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 773 } 774 775 for (i=0; i<A->m; i++) { 776 v = a->v + i; 777 for (j=0; j<A->n; j++) { 778 #if defined(PETSC_USE_COMPLEX) 779 if (allreal) { 780 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); 781 } else { 782 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 783 } 784 #else 785 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); 786 #endif 787 v += a->lda; 788 } 789 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 790 } 791 if (format == PETSC_VIEWER_ASCII_MATLAB) { 792 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 793 } 794 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 795 } 796 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 797 PetscFunctionReturn(0); 798 } 799 800 #undef __FUNCT__ 801 #define __FUNCT__ "MatView_SeqDense_Binary" 802 static int MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 803 { 804 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 805 int ict,j,n = A->n,m = A->m,i,fd,*col_lens,ierr,nz = m*n; 806 PetscScalar *v,*anonz,*vals; 807 PetscViewerFormat format; 808 809 PetscFunctionBegin; 810 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 811 812 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 813 if (format == PETSC_VIEWER_BINARY_NATIVE) { 814 /* store the matrix as a dense matrix */ 815 ierr = PetscMalloc(4*sizeof(int),&col_lens);CHKERRQ(ierr); 816 col_lens[0] = MAT_FILE_COOKIE; 817 col_lens[1] = m; 818 col_lens[2] = n; 819 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 820 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 821 ierr = PetscFree(col_lens);CHKERRQ(ierr); 822 823 /* write out matrix, by rows */ 824 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 825 v = a->v; 826 for (i=0; i<m; i++) { 827 for (j=0; j<n; j++) { 828 vals[i + j*m] = *v++; 829 } 830 } 831 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 832 ierr = PetscFree(vals);CHKERRQ(ierr); 833 } else { 834 ierr = PetscMalloc((4+nz)*sizeof(int),&col_lens);CHKERRQ(ierr); 835 col_lens[0] = MAT_FILE_COOKIE; 836 col_lens[1] = m; 837 col_lens[2] = n; 838 col_lens[3] = nz; 839 840 /* store lengths of each row and write (including header) to file */ 841 for (i=0; i<m; i++) col_lens[4+i] = n; 842 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 843 844 /* Possibly should write in smaller increments, not whole matrix at once? */ 845 /* store column indices (zero start index) */ 846 ict = 0; 847 for (i=0; i<m; i++) { 848 for (j=0; j<n; j++) col_lens[ict++] = j; 849 } 850 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 851 ierr = PetscFree(col_lens);CHKERRQ(ierr); 852 853 /* store nonzero values */ 854 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 855 ict = 0; 856 for (i=0; i<m; i++) { 857 v = a->v + i; 858 for (j=0; j<n; j++) { 859 anonz[ict++] = *v; v += a->lda; 860 } 861 } 862 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 863 ierr = PetscFree(anonz);CHKERRQ(ierr); 864 } 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNCT__ 869 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 870 int MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 871 { 872 Mat A = (Mat) Aa; 873 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 874 int m = A->m,n = A->n,color,i,j,ierr; 875 PetscScalar *v = a->v; 876 PetscViewer viewer; 877 PetscDraw popup; 878 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 879 PetscViewerFormat format; 880 881 PetscFunctionBegin; 882 883 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 884 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 885 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 886 887 /* Loop over matrix elements drawing boxes */ 888 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 889 /* Blue for negative and Red for positive */ 890 color = PETSC_DRAW_BLUE; 891 for(j = 0; j < n; j++) { 892 x_l = j; 893 x_r = x_l + 1.0; 894 for(i = 0; i < m; i++) { 895 y_l = m - i - 1.0; 896 y_r = y_l + 1.0; 897 #if defined(PETSC_USE_COMPLEX) 898 if (PetscRealPart(v[j*m+i]) > 0.) { 899 color = PETSC_DRAW_RED; 900 } else if (PetscRealPart(v[j*m+i]) < 0.) { 901 color = PETSC_DRAW_BLUE; 902 } else { 903 continue; 904 } 905 #else 906 if (v[j*m+i] > 0.) { 907 color = PETSC_DRAW_RED; 908 } else if (v[j*m+i] < 0.) { 909 color = PETSC_DRAW_BLUE; 910 } else { 911 continue; 912 } 913 #endif 914 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 915 } 916 } 917 } else { 918 /* use contour shading to indicate magnitude of values */ 919 /* first determine max of all nonzero values */ 920 for(i = 0; i < m*n; i++) { 921 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 922 } 923 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 924 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 925 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 926 for(j = 0; j < n; j++) { 927 x_l = j; 928 x_r = x_l + 1.0; 929 for(i = 0; i < m; i++) { 930 y_l = m - i - 1.0; 931 y_r = y_l + 1.0; 932 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 933 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 934 } 935 } 936 } 937 PetscFunctionReturn(0); 938 } 939 940 #undef __FUNCT__ 941 #define __FUNCT__ "MatView_SeqDense_Draw" 942 int MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 943 { 944 PetscDraw draw; 945 PetscTruth isnull; 946 PetscReal xr,yr,xl,yl,h,w; 947 int ierr; 948 949 PetscFunctionBegin; 950 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 951 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 952 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 953 954 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 955 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 956 xr += w; yr += h; xl = -w; yl = -h; 957 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 958 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 959 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 #undef __FUNCT__ 964 #define __FUNCT__ "MatView_SeqDense" 965 int MatView_SeqDense(Mat A,PetscViewer viewer) 966 { 967 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 968 int ierr; 969 PetscTruth issocket,isascii,isbinary,isdraw; 970 971 PetscFunctionBegin; 972 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 973 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 974 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 975 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 976 977 if (issocket) { 978 if (a->lda>A->m) SETERRQ(1,"Case can not handle LDA"); 979 ierr = PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);CHKERRQ(ierr); 980 } else if (isascii) { 981 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 982 } else if (isbinary) { 983 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 984 } else if (isdraw) { 985 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 986 } else { 987 SETERRQ1(1,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 988 } 989 PetscFunctionReturn(0); 990 } 991 992 #undef __FUNCT__ 993 #define __FUNCT__ "MatDestroy_SeqDense" 994 int MatDestroy_SeqDense(Mat mat) 995 { 996 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 997 int ierr; 998 999 PetscFunctionBegin; 1000 #if defined(PETSC_USE_LOG) 1001 PetscLogObjectState((PetscObject)mat,"Rows %d Cols %d",mat->m,mat->n); 1002 #endif 1003 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 1004 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1005 ierr = PetscFree(l);CHKERRQ(ierr); 1006 PetscFunctionReturn(0); 1007 } 1008 1009 #undef __FUNCT__ 1010 #define __FUNCT__ "MatTranspose_SeqDense" 1011 int MatTranspose_SeqDense(Mat A,Mat *matout) 1012 { 1013 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1014 int k,j,m,n,M,ierr; 1015 PetscScalar *v,tmp; 1016 1017 PetscFunctionBegin; 1018 v = mat->v; m = A->m; M = mat->lda; n = A->n; 1019 if (!matout) { /* in place transpose */ 1020 if (m != n) { 1021 SETERRQ(1,"Can not transpose non-square matrix in place"); 1022 } else { 1023 for (j=0; j<m; j++) { 1024 for (k=0; k<j; k++) { 1025 tmp = v[j + k*M]; 1026 v[j + k*M] = v[k + j*M]; 1027 v[k + j*M] = tmp; 1028 } 1029 } 1030 } 1031 } else { /* out-of-place transpose */ 1032 Mat tmat; 1033 Mat_SeqDense *tmatd; 1034 PetscScalar *v2; 1035 1036 ierr = MatCreateSeqDense(A->comm,A->n,A->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 1037 tmatd = (Mat_SeqDense*)tmat->data; 1038 v = mat->v; v2 = tmatd->v; 1039 for (j=0; j<n; j++) { 1040 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1041 } 1042 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1043 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1044 *matout = tmat; 1045 } 1046 PetscFunctionReturn(0); 1047 } 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "MatEqual_SeqDense" 1051 int MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1052 { 1053 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1054 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1055 int i,j; 1056 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1057 1058 PetscFunctionBegin; 1059 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1060 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1061 for (i=0; i<A1->m; i++) { 1062 v1 = mat1->v+i; v2 = mat2->v+i; 1063 for (j=0; j<A1->n; j++) { 1064 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1065 v1 += mat1->lda; v2 += mat2->lda; 1066 } 1067 } 1068 *flg = PETSC_TRUE; 1069 PetscFunctionReturn(0); 1070 } 1071 1072 #undef __FUNCT__ 1073 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1074 int MatGetDiagonal_SeqDense(Mat A,Vec v) 1075 { 1076 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1077 int ierr,i,n,len; 1078 PetscScalar *x,zero = 0.0; 1079 1080 PetscFunctionBegin; 1081 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1082 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1083 ierr = VecGetArrayFast(v,&x);CHKERRQ(ierr); 1084 len = PetscMin(A->m,A->n); 1085 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1086 for (i=0; i<len; i++) { 1087 x[i] = mat->v[i*mat->lda + i]; 1088 } 1089 ierr = VecRestoreArrayFast(v,&x);CHKERRQ(ierr); 1090 PetscFunctionReturn(0); 1091 } 1092 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1095 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1096 { 1097 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1098 PetscScalar *l,*r,x,*v; 1099 int ierr,i,j,m = A->m,n = A->n; 1100 1101 PetscFunctionBegin; 1102 if (ll) { 1103 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1104 ierr = VecGetArrayFast(ll,&l);CHKERRQ(ierr); 1105 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1106 for (i=0; i<m; i++) { 1107 x = l[i]; 1108 v = mat->v + i; 1109 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1110 } 1111 ierr = VecRestoreArrayFast(ll,&l);CHKERRQ(ierr); 1112 PetscLogFlops(n*m); 1113 } 1114 if (rr) { 1115 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1116 ierr = VecGetArrayFast(rr,&r);CHKERRQ(ierr); 1117 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1118 for (i=0; i<n; i++) { 1119 x = r[i]; 1120 v = mat->v + i*m; 1121 for (j=0; j<m; j++) { (*v++) *= x;} 1122 } 1123 ierr = VecRestoreArrayFast(rr,&r);CHKERRQ(ierr); 1124 PetscLogFlops(n*m); 1125 } 1126 PetscFunctionReturn(0); 1127 } 1128 1129 #undef __FUNCT__ 1130 #define __FUNCT__ "MatNorm_SeqDense" 1131 int MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1132 { 1133 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1134 PetscScalar *v = mat->v; 1135 PetscReal sum = 0.0; 1136 int lda=mat->lda,m=A->m,i,j; 1137 1138 PetscFunctionBegin; 1139 if (type == NORM_FROBENIUS) { 1140 if (lda>m) { 1141 for (j=0; j<A->n; j++) { 1142 v = mat->v+j*lda; 1143 for (i=0; i<m; i++) { 1144 #if defined(PETSC_USE_COMPLEX) 1145 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1146 #else 1147 sum += (*v)*(*v); v++; 1148 #endif 1149 } 1150 } 1151 } else { 1152 for (i=0; i<A->n*A->m; i++) { 1153 #if defined(PETSC_USE_COMPLEX) 1154 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1155 #else 1156 sum += (*v)*(*v); v++; 1157 #endif 1158 } 1159 } 1160 *nrm = sqrt(sum); 1161 PetscLogFlops(2*A->n*A->m); 1162 } else if (type == NORM_1) { 1163 *nrm = 0.0; 1164 for (j=0; j<A->n; j++) { 1165 v = mat->v + j*mat->lda; 1166 sum = 0.0; 1167 for (i=0; i<A->m; i++) { 1168 sum += PetscAbsScalar(*v); v++; 1169 } 1170 if (sum > *nrm) *nrm = sum; 1171 } 1172 PetscLogFlops(A->n*A->m); 1173 } else if (type == NORM_INFINITY) { 1174 *nrm = 0.0; 1175 for (j=0; j<A->m; j++) { 1176 v = mat->v + j; 1177 sum = 0.0; 1178 for (i=0; i<A->n; i++) { 1179 sum += PetscAbsScalar(*v); v += mat->lda; 1180 } 1181 if (sum > *nrm) *nrm = sum; 1182 } 1183 PetscLogFlops(A->n*A->m); 1184 } else { 1185 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1186 } 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNCT__ 1191 #define __FUNCT__ "MatSetOption_SeqDense" 1192 int MatSetOption_SeqDense(Mat A,MatOption op) 1193 { 1194 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1195 1196 PetscFunctionBegin; 1197 switch (op) { 1198 case MAT_ROW_ORIENTED: 1199 aij->roworiented = PETSC_TRUE; 1200 break; 1201 case MAT_COLUMN_ORIENTED: 1202 aij->roworiented = PETSC_FALSE; 1203 break; 1204 case MAT_ROWS_SORTED: 1205 case MAT_ROWS_UNSORTED: 1206 case MAT_COLUMNS_SORTED: 1207 case MAT_COLUMNS_UNSORTED: 1208 case MAT_NO_NEW_NONZERO_LOCATIONS: 1209 case MAT_YES_NEW_NONZERO_LOCATIONS: 1210 case MAT_NEW_NONZERO_LOCATION_ERR: 1211 case MAT_NO_NEW_DIAGONALS: 1212 case MAT_YES_NEW_DIAGONALS: 1213 case MAT_IGNORE_OFF_PROC_ENTRIES: 1214 case MAT_USE_HASH_TABLE: 1215 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1216 break; 1217 case MAT_SYMMETRIC: 1218 case MAT_STRUCTURALLY_SYMMETRIC: 1219 case MAT_NOT_SYMMETRIC: 1220 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 1221 case MAT_HERMITIAN: 1222 case MAT_NOT_HERMITIAN: 1223 case MAT_SYMMETRY_ETERNAL: 1224 case MAT_NOT_SYMMETRY_ETERNAL: 1225 break; 1226 default: 1227 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1228 } 1229 PetscFunctionReturn(0); 1230 } 1231 1232 #undef __FUNCT__ 1233 #define __FUNCT__ "MatZeroEntries_SeqDense" 1234 int MatZeroEntries_SeqDense(Mat A) 1235 { 1236 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1237 int lda=l->lda,m=A->m,j, ierr; 1238 1239 PetscFunctionBegin; 1240 if (lda>m) { 1241 for (j=0; j<A->n; j++) { 1242 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar)); CHKERRQ(ierr); 1243 } 1244 } else { 1245 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1246 } 1247 PetscFunctionReturn(0); 1248 } 1249 1250 #undef __FUNCT__ 1251 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1252 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1253 { 1254 PetscFunctionBegin; 1255 *bs = 1; 1256 PetscFunctionReturn(0); 1257 } 1258 1259 #undef __FUNCT__ 1260 #define __FUNCT__ "MatZeroRows_SeqDense" 1261 int MatZeroRows_SeqDense(Mat A,IS is,const PetscScalar *diag) 1262 { 1263 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1264 int n = A->n,i,j,ierr,N,*rows; 1265 PetscScalar *slot; 1266 1267 PetscFunctionBegin; 1268 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1269 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1270 for (i=0; i<N; i++) { 1271 slot = l->v + rows[i]; 1272 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1273 } 1274 if (diag) { 1275 for (i=0; i<N; i++) { 1276 slot = l->v + (n+1)*rows[i]; 1277 *slot = *diag; 1278 } 1279 } 1280 ISRestoreIndices(is,&rows); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "MatGetArray_SeqDense" 1286 int MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1287 { 1288 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1289 1290 PetscFunctionBegin; 1291 *array = mat->v; 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatRestoreArray_SeqDense" 1297 int MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1298 { 1299 PetscFunctionBegin; 1300 *array = 0; /* user cannot accidently use the array later */ 1301 PetscFunctionReturn(0); 1302 } 1303 1304 #undef __FUNCT__ 1305 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1306 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1307 { 1308 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1309 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1310 PetscScalar *av,*bv,*v = mat->v; 1311 Mat newmat; 1312 1313 PetscFunctionBegin; 1314 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1315 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1316 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1317 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1318 1319 /* Check submatrixcall */ 1320 if (scall == MAT_REUSE_MATRIX) { 1321 int n_cols,n_rows; 1322 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1323 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1324 newmat = *B; 1325 } else { 1326 /* Create and fill new matrix */ 1327 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1328 } 1329 1330 /* Now extract the data pointers and do the copy,column at a time */ 1331 bv = ((Mat_SeqDense*)newmat->data)->v; 1332 1333 for (i=0; i<ncols; i++) { 1334 av = v + m*icol[i]; 1335 for (j=0; j<nrows; j++) { 1336 *bv++ = av[irow[j]]; 1337 } 1338 } 1339 1340 /* Assemble the matrices so that the correct flags are set */ 1341 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1342 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1343 1344 /* Free work space */ 1345 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1346 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1347 *B = newmat; 1348 PetscFunctionReturn(0); 1349 } 1350 1351 #undef __FUNCT__ 1352 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1353 int MatGetSubMatrices_SeqDense(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1354 { 1355 int ierr,i; 1356 1357 PetscFunctionBegin; 1358 if (scall == MAT_INITIAL_MATRIX) { 1359 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1360 } 1361 1362 for (i=0; i<n; i++) { 1363 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1364 } 1365 PetscFunctionReturn(0); 1366 } 1367 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "MatCopy_SeqDense" 1370 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1371 { 1372 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1373 int lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j,ierr; 1374 1375 PetscFunctionBegin; 1376 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1377 if (A->ops->copy != B->ops->copy) { 1378 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1379 PetscFunctionReturn(0); 1380 } 1381 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1382 if (lda1>m || lda2>m) { 1383 for (j=0; j<n; j++) { 1384 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1385 } 1386 } else { 1387 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1388 } 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1394 int MatSetUpPreallocation_SeqDense(Mat A) 1395 { 1396 int ierr; 1397 1398 PetscFunctionBegin; 1399 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1400 PetscFunctionReturn(0); 1401 } 1402 1403 /* -------------------------------------------------------------------*/ 1404 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1405 MatGetRow_SeqDense, 1406 MatRestoreRow_SeqDense, 1407 MatMult_SeqDense, 1408 /* 4*/ MatMultAdd_SeqDense, 1409 MatMultTranspose_SeqDense, 1410 MatMultTransposeAdd_SeqDense, 1411 MatSolve_SeqDense, 1412 MatSolveAdd_SeqDense, 1413 MatSolveTranspose_SeqDense, 1414 /*10*/ MatSolveTransposeAdd_SeqDense, 1415 MatLUFactor_SeqDense, 1416 MatCholeskyFactor_SeqDense, 1417 MatRelax_SeqDense, 1418 MatTranspose_SeqDense, 1419 /*15*/ MatGetInfo_SeqDense, 1420 MatEqual_SeqDense, 1421 MatGetDiagonal_SeqDense, 1422 MatDiagonalScale_SeqDense, 1423 MatNorm_SeqDense, 1424 /*20*/ 0, 1425 0, 1426 0, 1427 MatSetOption_SeqDense, 1428 MatZeroEntries_SeqDense, 1429 /*25*/ MatZeroRows_SeqDense, 1430 MatLUFactorSymbolic_SeqDense, 1431 MatLUFactorNumeric_SeqDense, 1432 MatCholeskyFactorSymbolic_SeqDense, 1433 MatCholeskyFactorNumeric_SeqDense, 1434 /*30*/ MatSetUpPreallocation_SeqDense, 1435 0, 1436 0, 1437 MatGetArray_SeqDense, 1438 MatRestoreArray_SeqDense, 1439 /*35*/ MatDuplicate_SeqDense, 1440 0, 1441 0, 1442 0, 1443 0, 1444 /*40*/ MatAXPY_SeqDense, 1445 MatGetSubMatrices_SeqDense, 1446 0, 1447 MatGetValues_SeqDense, 1448 MatCopy_SeqDense, 1449 /*45*/ 0, 1450 MatScale_SeqDense, 1451 0, 1452 0, 1453 0, 1454 /*50*/ MatGetBlockSize_SeqDense, 1455 0, 1456 0, 1457 0, 1458 0, 1459 /*55*/ 0, 1460 0, 1461 0, 1462 0, 1463 0, 1464 /*60*/ 0, 1465 MatDestroy_SeqDense, 1466 MatView_SeqDense, 1467 MatGetPetscMaps_Petsc, 1468 0, 1469 /*65*/ 0, 1470 0, 1471 0, 1472 0, 1473 0, 1474 /*70*/ 0, 1475 0, 1476 0, 1477 0, 1478 0, 1479 /*75*/ 0, 1480 0, 1481 0, 1482 0, 1483 0, 1484 /*80*/ 0, 1485 0, 1486 0, 1487 0, 1488 /*85*/ MatLoad_SeqDense}; 1489 1490 #undef __FUNCT__ 1491 #define __FUNCT__ "MatCreateSeqDense" 1492 /*@C 1493 MatCreateSeqDense - Creates a sequential dense matrix that 1494 is stored in column major order (the usual Fortran 77 manner). Many 1495 of the matrix operations use the BLAS and LAPACK routines. 1496 1497 Collective on MPI_Comm 1498 1499 Input Parameters: 1500 + comm - MPI communicator, set to PETSC_COMM_SELF 1501 . m - number of rows 1502 . n - number of columns 1503 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1504 to control all matrix memory allocation. 1505 1506 Output Parameter: 1507 . A - the matrix 1508 1509 Notes: 1510 The data input variable is intended primarily for Fortran programmers 1511 who wish to allocate their own matrix memory space. Most users should 1512 set data=PETSC_NULL. 1513 1514 Level: intermediate 1515 1516 .keywords: dense, matrix, LAPACK, BLAS 1517 1518 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1519 @*/ 1520 int MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1521 { 1522 int ierr; 1523 1524 PetscFunctionBegin; 1525 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1526 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1527 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1528 PetscFunctionReturn(0); 1529 } 1530 1531 #undef __FUNCT__ 1532 #define __FUNCT__ "MatSeqDensePreallocation" 1533 /*@C 1534 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1535 1536 Collective on MPI_Comm 1537 1538 Input Parameters: 1539 + A - the matrix 1540 - data - the array (or PETSC_NULL) 1541 1542 Notes: 1543 The data input variable is intended primarily for Fortran programmers 1544 who wish to allocate their own matrix memory space. Most users should 1545 set data=PETSC_NULL. 1546 1547 Level: intermediate 1548 1549 .keywords: dense, matrix, LAPACK, BLAS 1550 1551 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1552 @*/ 1553 int MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1554 { 1555 int ierr,(*f)(Mat,PetscScalar[]); 1556 1557 PetscFunctionBegin; 1558 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1559 if (f) { 1560 ierr = (*f)(B,data);CHKERRQ(ierr); 1561 } 1562 PetscFunctionReturn(0); 1563 } 1564 1565 EXTERN_C_BEGIN 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1568 int MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1569 { 1570 Mat_SeqDense *b; 1571 int ierr; 1572 1573 PetscFunctionBegin; 1574 B->preallocated = PETSC_TRUE; 1575 b = (Mat_SeqDense*)B->data; 1576 if (!data) { 1577 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1578 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1579 b->user_alloc = PETSC_FALSE; 1580 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1581 } else { /* user-allocated storage */ 1582 b->v = data; 1583 b->user_alloc = PETSC_TRUE; 1584 } 1585 PetscFunctionReturn(0); 1586 } 1587 EXTERN_C_END 1588 1589 #undef __FUNCT__ 1590 #define __FUNCT__ "MatSeqDenseSetLDA" 1591 /*@C 1592 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1593 1594 Input parameter: 1595 + A - the matrix 1596 - lda - the leading dimension 1597 1598 Notes: 1599 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1600 it asserts that the preallocation has a leading dimension (the LDA parameter 1601 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1602 1603 Level: intermediate 1604 1605 .keywords: dense, matrix, LAPACK, BLAS 1606 1607 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1608 @*/ 1609 int MatSeqDenseSetLDA(Mat B,int lda) 1610 { 1611 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1612 PetscFunctionBegin; 1613 if (lda<B->m) SETERRQ(1,"LDA must be at least matrix i dimension"); 1614 b->lda = lda; 1615 PetscFunctionReturn(0); 1616 } 1617 1618 /*MC 1619 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1620 1621 Options Database Keys: 1622 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1623 1624 Level: beginner 1625 1626 .seealso: MatCreateSeqDense 1627 M*/ 1628 1629 EXTERN_C_BEGIN 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatCreate_SeqDense" 1632 int MatCreate_SeqDense(Mat B) 1633 { 1634 Mat_SeqDense *b; 1635 int ierr,size; 1636 1637 PetscFunctionBegin; 1638 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1639 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1640 1641 B->m = B->M = PetscMax(B->m,B->M); 1642 B->n = B->N = PetscMax(B->n,B->N); 1643 1644 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1645 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1646 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1647 B->factor = 0; 1648 B->mapping = 0; 1649 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1650 B->data = (void*)b; 1651 1652 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1653 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1654 1655 b->pivots = 0; 1656 b->roworiented = PETSC_TRUE; 1657 b->v = 0; 1658 b->lda = B->m; 1659 1660 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1661 "MatSeqDenseSetPreallocation_SeqDense", 1662 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1663 PetscFunctionReturn(0); 1664 } 1665 EXTERN_C_END 1666