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