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