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