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