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