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