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