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