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