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