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