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_DETAIL) { 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 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1180 break; 1181 default: 1182 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1183 } 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatZeroEntries_SeqDense" 1189 int MatZeroEntries_SeqDense(Mat A) 1190 { 1191 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1192 int ierr; 1193 1194 PetscFunctionBegin; 1195 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1201 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1202 { 1203 PetscFunctionBegin; 1204 *bs = 1; 1205 PetscFunctionReturn(0); 1206 } 1207 1208 #undef __FUNCT__ 1209 #define __FUNCT__ "MatZeroRows_SeqDense" 1210 int MatZeroRows_SeqDense(Mat A,IS is,PetscScalar *diag) 1211 { 1212 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1213 int n = A->n,i,j,ierr,N,*rows; 1214 PetscScalar *slot; 1215 1216 PetscFunctionBegin; 1217 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1218 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1219 for (i=0; i<N; i++) { 1220 slot = l->v + rows[i]; 1221 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1222 } 1223 if (diag) { 1224 for (i=0; i<N; i++) { 1225 slot = l->v + (n+1)*rows[i]; 1226 *slot = *diag; 1227 } 1228 } 1229 ISRestoreIndices(is,&rows); 1230 PetscFunctionReturn(0); 1231 } 1232 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "MatGetArray_SeqDense" 1235 int MatGetArray_SeqDense(Mat A,PetscScalar **array) 1236 { 1237 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1238 1239 PetscFunctionBegin; 1240 *array = mat->v; 1241 PetscFunctionReturn(0); 1242 } 1243 1244 #undef __FUNCT__ 1245 #define __FUNCT__ "MatRestoreArray_SeqDense" 1246 int MatRestoreArray_SeqDense(Mat A,PetscScalar **array) 1247 { 1248 PetscFunctionBegin; 1249 *array = 0; /* user cannot accidently use the array later */ 1250 PetscFunctionReturn(0); 1251 } 1252 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1255 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1256 { 1257 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1258 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1259 PetscScalar *av,*bv,*v = mat->v; 1260 Mat newmat; 1261 1262 PetscFunctionBegin; 1263 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1264 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1265 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1266 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1267 1268 /* Check submatrixcall */ 1269 if (scall == MAT_REUSE_MATRIX) { 1270 int n_cols,n_rows; 1271 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1272 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1273 newmat = *B; 1274 } else { 1275 /* Create and fill new matrix */ 1276 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1277 } 1278 1279 /* Now extract the data pointers and do the copy,column at a time */ 1280 bv = ((Mat_SeqDense*)newmat->data)->v; 1281 1282 for (i=0; i<ncols; i++) { 1283 av = v + m*icol[i]; 1284 for (j=0; j<nrows; j++) { 1285 *bv++ = av[irow[j]]; 1286 } 1287 } 1288 1289 /* Assemble the matrices so that the correct flags are set */ 1290 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1291 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1292 1293 /* Free work space */ 1294 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1295 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1296 *B = newmat; 1297 PetscFunctionReturn(0); 1298 } 1299 1300 #undef __FUNCT__ 1301 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1302 int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1303 { 1304 int ierr,i; 1305 1306 PetscFunctionBegin; 1307 if (scall == MAT_INITIAL_MATRIX) { 1308 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1309 } 1310 1311 for (i=0; i<n; i++) { 1312 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1313 } 1314 PetscFunctionReturn(0); 1315 } 1316 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "MatCopy_SeqDense" 1319 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1320 { 1321 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1322 int ierr; 1323 PetscTruth flag; 1324 1325 PetscFunctionBegin; 1326 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flag);CHKERRQ(ierr); 1327 if (!flag) { 1328 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 if (A->m != B->m || A->n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1332 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1333 PetscFunctionReturn(0); 1334 } 1335 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1338 int MatSetUpPreallocation_SeqDense(Mat A) 1339 { 1340 int ierr; 1341 1342 PetscFunctionBegin; 1343 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1344 PetscFunctionReturn(0); 1345 } 1346 1347 /* -------------------------------------------------------------------*/ 1348 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1349 MatGetRow_SeqDense, 1350 MatRestoreRow_SeqDense, 1351 MatMult_SeqDense, 1352 MatMultAdd_SeqDense, 1353 MatMultTranspose_SeqDense, 1354 MatMultTransposeAdd_SeqDense, 1355 MatSolve_SeqDense, 1356 MatSolveAdd_SeqDense, 1357 MatSolveTranspose_SeqDense, 1358 MatSolveTransposeAdd_SeqDense, 1359 MatLUFactor_SeqDense, 1360 MatCholeskyFactor_SeqDense, 1361 MatRelax_SeqDense, 1362 MatTranspose_SeqDense, 1363 MatGetInfo_SeqDense, 1364 MatEqual_SeqDense, 1365 MatGetDiagonal_SeqDense, 1366 MatDiagonalScale_SeqDense, 1367 MatNorm_SeqDense, 1368 0, 1369 0, 1370 0, 1371 MatSetOption_SeqDense, 1372 MatZeroEntries_SeqDense, 1373 MatZeroRows_SeqDense, 1374 MatLUFactorSymbolic_SeqDense, 1375 MatLUFactorNumeric_SeqDense, 1376 MatCholeskyFactorSymbolic_SeqDense, 1377 MatCholeskyFactorNumeric_SeqDense, 1378 MatSetUpPreallocation_SeqDense, 1379 0, 1380 0, 1381 MatGetArray_SeqDense, 1382 MatRestoreArray_SeqDense, 1383 MatDuplicate_SeqDense, 1384 0, 1385 0, 1386 0, 1387 0, 1388 MatAXPY_SeqDense, 1389 MatGetSubMatrices_SeqDense, 1390 0, 1391 MatGetValues_SeqDense, 1392 MatCopy_SeqDense, 1393 0, 1394 MatScale_SeqDense, 1395 0, 1396 0, 1397 0, 1398 MatGetBlockSize_SeqDense, 1399 0, 1400 0, 1401 0, 1402 0, 1403 0, 1404 0, 1405 0, 1406 0, 1407 0, 1408 0, 1409 MatDestroy_SeqDense, 1410 MatView_SeqDense, 1411 MatGetPetscMaps_Petsc}; 1412 1413 #undef __FUNCT__ 1414 #define __FUNCT__ "MatCreateSeqDense" 1415 /*@C 1416 MatCreateSeqDense - Creates a sequential dense matrix that 1417 is stored in column major order (the usual Fortran 77 manner). Many 1418 of the matrix operations use the BLAS and LAPACK routines. 1419 1420 Collective on MPI_Comm 1421 1422 Input Parameters: 1423 + comm - MPI communicator, set to PETSC_COMM_SELF 1424 . m - number of rows 1425 . n - number of columns 1426 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1427 to control all matrix memory allocation. 1428 1429 Output Parameter: 1430 . A - the matrix 1431 1432 Notes: 1433 The data input variable is intended primarily for Fortran programmers 1434 who wish to allocate their own matrix memory space. Most users should 1435 set data=PETSC_NULL. 1436 1437 Level: intermediate 1438 1439 .keywords: dense, matrix, LAPACK, BLAS 1440 1441 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1442 @*/ 1443 int MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1444 { 1445 int ierr; 1446 1447 PetscFunctionBegin; 1448 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1449 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1450 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 #undef __FUNCT__ 1455 #define __FUNCT__ "MatSeqDensePreallocation" 1456 /*@C 1457 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1458 1459 Collective on MPI_Comm 1460 1461 Input Parameters: 1462 + A - the matrix 1463 - data - the array (or PETSC_NULL) 1464 1465 Notes: 1466 The data input variable is intended primarily for Fortran programmers 1467 who wish to allocate their own matrix memory space. Most users should 1468 set data=PETSC_NULL. 1469 1470 Level: intermediate 1471 1472 .keywords: dense, matrix, LAPACK, BLAS 1473 1474 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1475 @*/ 1476 int MatSeqDenseSetPreallocation(Mat B,PetscScalar *data) 1477 { 1478 Mat_SeqDense *b; 1479 int ierr; 1480 PetscTruth flg2; 1481 1482 PetscFunctionBegin; 1483 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flg2);CHKERRQ(ierr); 1484 if (!flg2) PetscFunctionReturn(0); 1485 B->preallocated = PETSC_TRUE; 1486 b = (Mat_SeqDense*)B->data; 1487 if (!data) { 1488 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1489 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1490 b->user_alloc = PETSC_FALSE; 1491 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1492 } else { /* user-allocated storage */ 1493 b->v = data; 1494 b->user_alloc = PETSC_TRUE; 1495 } 1496 PetscFunctionReturn(0); 1497 } 1498 1499 EXTERN_C_BEGIN 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "MatCreate_SeqDense" 1502 int MatCreate_SeqDense(Mat B) 1503 { 1504 Mat_SeqDense *b; 1505 int ierr,size; 1506 1507 PetscFunctionBegin; 1508 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1509 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1510 1511 B->m = B->M = PetscMax(B->m,B->M); 1512 B->n = B->N = PetscMax(B->n,B->N); 1513 1514 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1515 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1516 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1517 B->factor = 0; 1518 B->mapping = 0; 1519 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1520 B->data = (void*)b; 1521 1522 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1523 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1524 1525 b->pivots = 0; 1526 b->roworiented = PETSC_TRUE; 1527 b->v = 0; 1528 1529 PetscFunctionReturn(0); 1530 } 1531 EXTERN_C_END 1532 1533 1534 1535 1536