1 /*$Id: dense.c,v 1.190 2000/10/24 20:25:29 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__ /*<a name=""></a>*/"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 PLogFlops(2*N-1); 19 PetscFunctionReturn(0); 20 } 21 22 #undef __FUNC__ 23 #define __FUNC__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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 PLogFlops(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__ /*<a name=""></a>*/"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; 74 75 PetscFunctionBegin; 76 if (!mat->pivots) { 77 mat->pivots = (int*)PetscMalloc((A->m+1)*sizeof(int));CHKPTRQ(mat->pivots); 78 PLogObjectMemory(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 PLogFlops((2*A->n*A->n*A->n)/3); 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNC__ 90 #define __FUNC__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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 PLogObjectMemory(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 PLogFlops((A->n*A->n*A->n)/3); 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNC__ 178 #define __FUNC__ /*<a name=""></a>*/"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 PLogFlops(2*A->n*A->n - A->n); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNC__ 204 #define __FUNC__ /*<a name=""></a>*/"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 PLogFlops(2*A->n*A->n - A->n); 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNC__ 230 #define __FUNC__ /*<a name=""></a>*/"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 PLogObjectParent(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 PLogFlops(2*A->n*A->n); 260 PetscFunctionReturn(0); 261 } 262 263 #undef __FUNC__ 264 #define __FUNC__ /*<a name=""></a>*/"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 PLogObjectParent(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 PLogFlops(2*A->n*A->n); 298 PetscFunctionReturn(0); 299 } 300 /* ------------------------------------------------------------------*/ 301 #undef __FUNC__ 302 #define __FUNC__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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 PLogFlops(2*A->m*A->n - A->n); 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNC__ 382 #define __FUNC__ /*<a name=""></a>*/"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 PLogFlops(2*A->m*A->n - A->m); 397 PetscFunctionReturn(0); 398 } 399 400 #undef __FUNC__ 401 #define __FUNC__ /*<a name=""></a>*/"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 PLogFlops(2*A->m*A->n); 417 PetscFunctionReturn(0); 418 } 419 420 #undef __FUNC__ 421 #define __FUNC__ /*<a name=""></a>*/"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 PLogFlops(2*A->m*A->n); 438 PetscFunctionReturn(0); 439 } 440 441 /* -----------------------------------------------------------------*/ 442 #undef __FUNC__ 443 #define __FUNC__ /*<a name=""></a>*/"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; 449 450 PetscFunctionBegin; 451 *ncols = A->n; 452 if (cols) { 453 *cols = (int*)PetscMalloc((A->n+1)*sizeof(int));CHKPTRQ(*cols); 454 for (i=0; i<A->n; i++) (*cols)[i] = i; 455 } 456 if (vals) { 457 *vals = (Scalar*)PetscMalloc((A->n+1)*sizeof(Scalar));CHKPTRQ(*vals); 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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"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__ /*<a name=""></a>*/"MatLoad_SeqDense" 573 int MatLoad_SeqDense(Viewer 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 = ViewerBinaryGetDescriptor(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 w = (Scalar*)PetscMalloc((M*N+1)*sizeof(Scalar));CHKPTRQ(w); 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 rowlengths = (int*)PetscMalloc((M+1)*sizeof(int));CHKPTRQ(rowlengths); 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 cols = scols = (int*)PetscMalloc((nz+1)*sizeof(int));CHKPTRQ(cols); 622 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 623 vals = svals = (Scalar*)PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(vals); 624 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 625 626 /* insert into matrix */ 627 for (i=0; i<M; i++) { 628 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 629 svals += rowlengths[i]; scols += rowlengths[i]; 630 } 631 ierr = PetscFree(vals);CHKERRQ(ierr); 632 ierr = PetscFree(cols);CHKERRQ(ierr); 633 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 634 635 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 636 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 637 } 638 PetscFunctionReturn(0); 639 } 640 EXTERN_C_END 641 642 #include "petscsys.h" 643 644 #undef __FUNC__ 645 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_ASCII" 646 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 647 { 648 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 649 int ierr,i,j,format; 650 char *outputname; 651 Scalar *v; 652 653 PetscFunctionBegin; 654 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 655 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 656 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 657 PetscFunctionReturn(0); /* do nothing for now */ 658 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 659 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 660 for (i=0; i<A->m; i++) { 661 v = a->v + i; 662 ierr = ViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 663 for (j=0; j<A->n; j++) { 664 #if defined(PETSC_USE_COMPLEX) 665 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 666 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 667 } else if (PetscRealPart(*v)) { 668 ierr = ViewerASCIIPrintf(viewer," %d %g ",j,PetscRealPart(*v));CHKERRQ(ierr); 669 } 670 #else 671 if (*v) { 672 ierr = ViewerASCIIPrintf(viewer," %d %g ",j,*v);CHKERRQ(ierr); 673 } 674 #endif 675 v += A->m; 676 } 677 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 678 } 679 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 680 } else { 681 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 682 #if defined(PETSC_USE_COMPLEX) 683 PetscTruth allreal = PETSC_TRUE; 684 /* determine if matrix has all real values */ 685 v = a->v; 686 for (i=0; i<A->m*A->n; i++) { 687 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 688 } 689 #endif 690 if (format == VIEWER_FORMAT_ASCII_MATLAB) { 691 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 692 if (!outputname) outputname = "A"; 693 ierr = ViewerASCIIPrintf(viewer,"%% Size = %d %d \n",A->m,A->n);CHKERRQ(ierr); 694 ierr = ViewerASCIIPrintf(viewer,"%s = zeros(%d,%d);\n",outputname,A->m,A->n);CHKERRQ(ierr); 695 ierr = ViewerASCIIPrintf(viewer,"%s = [\n",outputname);CHKERRQ(ierr); 696 } 697 698 for (i=0; i<A->m; i++) { 699 v = a->v + i; 700 for (j=0; j<A->n; j++) { 701 #if defined(PETSC_USE_COMPLEX) 702 if (allreal) { 703 ierr = ViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); v += A->m; 704 } else { 705 ierr = ViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); v += A->m; 706 } 707 #else 708 ierr = ViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); v += A->m; 709 #endif 710 } 711 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 712 } 713 if (format == VIEWER_FORMAT_ASCII_MATLAB) { 714 ierr = ViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 715 } 716 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 717 } 718 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 719 PetscFunctionReturn(0); 720 } 721 722 #undef __FUNC__ 723 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Binary" 724 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 725 { 726 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 727 int ict,j,n = A->n,m = A->m,i,fd,*col_lens,ierr,nz = m*n; 728 int format; 729 Scalar *v,*anonz,*vals; 730 731 PetscFunctionBegin; 732 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 733 734 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 735 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 736 /* store the matrix as a dense matrix */ 737 col_lens = (int*)PetscMalloc(4*sizeof(int));CHKPTRQ(col_lens); 738 col_lens[0] = MAT_COOKIE; 739 col_lens[1] = m; 740 col_lens[2] = n; 741 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 742 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 743 ierr = PetscFree(col_lens);CHKERRQ(ierr); 744 745 /* write out matrix, by rows */ 746 vals = (Scalar*)PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(vals); 747 v = a->v; 748 for (i=0; i<m; i++) { 749 for (j=0; j<n; j++) { 750 vals[i + j*m] = *v++; 751 } 752 } 753 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 754 ierr = PetscFree(vals);CHKERRQ(ierr); 755 } else { 756 col_lens = (int*)PetscMalloc((4+nz)*sizeof(int));CHKPTRQ(col_lens); 757 col_lens[0] = MAT_COOKIE; 758 col_lens[1] = m; 759 col_lens[2] = n; 760 col_lens[3] = nz; 761 762 /* store lengths of each row and write (including header) to file */ 763 for (i=0; i<m; i++) col_lens[4+i] = n; 764 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 765 766 /* Possibly should write in smaller increments, not whole matrix at once? */ 767 /* store column indices (zero start index) */ 768 ict = 0; 769 for (i=0; i<m; i++) { 770 for (j=0; j<n; j++) col_lens[ict++] = j; 771 } 772 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 773 ierr = PetscFree(col_lens);CHKERRQ(ierr); 774 775 /* store nonzero values */ 776 anonz = (Scalar*)PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(anonz); 777 ict = 0; 778 for (i=0; i<m; i++) { 779 v = a->v + i; 780 for (j=0; j<n; j++) { 781 anonz[ict++] = *v; v += A->m; 782 } 783 } 784 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 785 ierr = PetscFree(anonz);CHKERRQ(ierr); 786 } 787 PetscFunctionReturn(0); 788 } 789 790 #undef __FUNC__ 791 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Draw_Zoom" 792 int MatView_SeqDense_Draw_Zoom(Draw draw,void *Aa) 793 { 794 Mat A = (Mat) Aa; 795 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 796 int m = A->m,n = A->n,format,color,i,j,ierr; 797 Scalar *v = a->v; 798 Viewer viewer; 799 Draw popup; 800 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 801 802 PetscFunctionBegin; 803 804 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 805 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 806 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 807 808 /* Loop over matrix elements drawing boxes */ 809 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 810 /* Blue for negative and Red for positive */ 811 color = DRAW_BLUE; 812 for(j = 0; j < n; j++) { 813 x_l = j; 814 x_r = x_l + 1.0; 815 for(i = 0; i < m; i++) { 816 y_l = m - i - 1.0; 817 y_r = y_l + 1.0; 818 #if defined(PETSC_USE_COMPLEX) 819 if (PetscRealPart(v[j*m+i]) > 0.) { 820 color = DRAW_RED; 821 } else if (PetscRealPart(v[j*m+i]) < 0.) { 822 color = DRAW_BLUE; 823 } else { 824 continue; 825 } 826 #else 827 if (v[j*m+i] > 0.) { 828 color = DRAW_RED; 829 } else if (v[j*m+i] < 0.) { 830 color = DRAW_BLUE; 831 } else { 832 continue; 833 } 834 #endif 835 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 836 } 837 } 838 } else { 839 /* use contour shading to indicate magnitude of values */ 840 /* first determine max of all nonzero values */ 841 for(i = 0; i < m*n; i++) { 842 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 843 } 844 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 845 ierr = DrawGetPopup(draw,&popup);CHKERRQ(ierr); 846 if (popup) {ierr = DrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 847 for(j = 0; j < n; j++) { 848 x_l = j; 849 x_r = x_l + 1.0; 850 for(i = 0; i < m; i++) { 851 y_l = m - i - 1.0; 852 y_r = y_l + 1.0; 853 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 854 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 855 } 856 } 857 } 858 PetscFunctionReturn(0); 859 } 860 861 #undef __FUNC__ 862 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Draw" 863 int MatView_SeqDense_Draw(Mat A,Viewer viewer) 864 { 865 Draw draw; 866 PetscTruth isnull; 867 PetscReal xr,yr,xl,yl,h,w; 868 int ierr; 869 870 PetscFunctionBegin; 871 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 872 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); 873 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 874 875 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 876 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 877 xr += w; yr += h; xl = -w; yl = -h; 878 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 879 ierr = DrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 880 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 881 PetscFunctionReturn(0); 882 } 883 884 #undef __FUNC__ 885 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense" 886 int MatView_SeqDense(Mat A,Viewer viewer) 887 { 888 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 889 int ierr; 890 PetscTruth issocket,isascii,isbinary,isdraw; 891 892 PetscFunctionBegin; 893 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 894 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 895 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 896 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 897 898 if (issocket) { 899 ierr = ViewerSocketPutScalar(viewer,A->m,A->n,a->v);CHKERRQ(ierr); 900 } else if (isascii) { 901 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 902 } else if (isbinary) { 903 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 904 } else if (isdraw) { 905 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 906 } else { 907 SETERRQ1(1,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 908 } 909 PetscFunctionReturn(0); 910 } 911 912 #undef __FUNC__ 913 #define __FUNC__ /*<a name=""></a>*/"MatDestroy_SeqDense" 914 int MatDestroy_SeqDense(Mat mat) 915 { 916 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 917 int ierr; 918 919 PetscFunctionBegin; 920 #if defined(PETSC_USE_LOG) 921 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",mat->m,mat->n); 922 #endif 923 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 924 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 925 ierr = PetscFree(l);CHKERRQ(ierr); 926 PetscFunctionReturn(0); 927 } 928 929 #undef __FUNC__ 930 #define __FUNC__ /*<a name=""></a>*/"MatTranspose_SeqDense" 931 int MatTranspose_SeqDense(Mat A,Mat *matout) 932 { 933 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 934 int k,j,m,n,ierr; 935 Scalar *v,tmp; 936 937 PetscFunctionBegin; 938 v = mat->v; m = A->m; n = A->n; 939 if (!matout) { /* in place transpose */ 940 if (m != n) { /* malloc temp to hold transpose */ 941 Scalar *w = (Scalar*)PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 942 for (j=0; j<m; j++) { 943 for (k=0; k<n; k++) { 944 w[k + j*n] = v[j + k*m]; 945 } 946 } 947 ierr = PetscMemcpy(v,w,m*n*sizeof(Scalar));CHKERRQ(ierr); 948 ierr = PetscFree(w);CHKERRQ(ierr); 949 } else { 950 for (j=0; j<m; j++) { 951 for (k=0; k<j; k++) { 952 tmp = v[j + k*n]; 953 v[j + k*n] = v[k + j*n]; 954 v[k + j*n] = tmp; 955 } 956 } 957 } 958 } else { /* out-of-place transpose */ 959 Mat tmat; 960 Mat_SeqDense *tmatd; 961 Scalar *v2; 962 ierr = MatCreateSeqDense(A->comm,A->n,A->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 963 tmatd = (Mat_SeqDense*)tmat->data; 964 v = mat->v; v2 = tmatd->v; 965 for (j=0; j<n; j++) { 966 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*m]; 967 } 968 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 969 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 970 *matout = tmat; 971 } 972 PetscFunctionReturn(0); 973 } 974 975 #undef __FUNC__ 976 #define __FUNC__ /*<a name=""></a>*/"MatEqual_SeqDense" 977 int MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 978 { 979 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 980 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 981 int ierr,i; 982 Scalar *v1 = mat1->v,*v2 = mat2->v; 983 PetscTruth flag; 984 985 PetscFunctionBegin; 986 ierr = PetscTypeCompare((PetscObject)A2,MATSEQDENSE,&flag);CHKERRQ(ierr); 987 if (!flag) SETERRQ(PETSC_ERR_SUP,"Matrices should be of type MATSEQDENSE"); 988 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 989 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 990 for (i=0; i<A1->m*A1->n; i++) { 991 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 992 v1++; v2++; 993 } 994 *flg = PETSC_TRUE; 995 PetscFunctionReturn(0); 996 } 997 998 #undef __FUNC__ 999 #define __FUNC__ /*<a name=""></a>*/"MatGetDiagonal_SeqDense" 1000 int MatGetDiagonal_SeqDense(Mat A,Vec v) 1001 { 1002 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1003 int ierr,i,n,len; 1004 Scalar *x,zero = 0.0; 1005 1006 PetscFunctionBegin; 1007 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1008 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1009 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1010 len = PetscMin(A->m,A->n); 1011 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1012 for (i=0; i<len; i++) { 1013 x[i] = mat->v[i*A->m + i]; 1014 } 1015 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #undef __FUNC__ 1020 #define __FUNC__ /*<a name=""></a>*/"MatDiagonalScale_SeqDense" 1021 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1022 { 1023 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1024 Scalar *l,*r,x,*v; 1025 int ierr,i,j,m = A->m,n = A->n; 1026 1027 PetscFunctionBegin; 1028 if (ll) { 1029 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1030 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1031 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1032 for (i=0; i<m; i++) { 1033 x = l[i]; 1034 v = mat->v + i; 1035 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1036 } 1037 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1038 PLogFlops(n*m); 1039 } 1040 if (rr) { 1041 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1042 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1043 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1044 for (i=0; i<n; i++) { 1045 x = r[i]; 1046 v = mat->v + i*m; 1047 for (j=0; j<m; j++) { (*v++) *= x;} 1048 } 1049 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1050 PLogFlops(n*m); 1051 } 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNC__ 1056 #define __FUNC__ /*<a name=""></a>*/"MatNorm_SeqDense" 1057 int MatNorm_SeqDense(Mat A,NormType type,PetscReal *norm) 1058 { 1059 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1060 Scalar *v = mat->v; 1061 PetscReal sum = 0.0; 1062 int i,j; 1063 1064 PetscFunctionBegin; 1065 if (type == NORM_FROBENIUS) { 1066 for (i=0; i<A->n*A->m; i++) { 1067 #if defined(PETSC_USE_COMPLEX) 1068 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1069 #else 1070 sum += (*v)*(*v); v++; 1071 #endif 1072 } 1073 *norm = sqrt(sum); 1074 PLogFlops(2*A->n*A->m); 1075 } else if (type == NORM_1) { 1076 *norm = 0.0; 1077 for (j=0; j<A->n; j++) { 1078 sum = 0.0; 1079 for (i=0; i<A->m; i++) { 1080 sum += PetscAbsScalar(*v); v++; 1081 } 1082 if (sum > *norm) *norm = sum; 1083 } 1084 PLogFlops(A->n*A->m); 1085 } else if (type == NORM_INFINITY) { 1086 *norm = 0.0; 1087 for (j=0; j<A->m; j++) { 1088 v = mat->v + j; 1089 sum = 0.0; 1090 for (i=0; i<A->n; i++) { 1091 sum += PetscAbsScalar(*v); v += A->m; 1092 } 1093 if (sum > *norm) *norm = sum; 1094 } 1095 PLogFlops(A->n*A->m); 1096 } else { 1097 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNC__ 1103 #define __FUNC__ /*<a name=""></a>*/"MatSetOption_SeqDense" 1104 int MatSetOption_SeqDense(Mat A,MatOption op) 1105 { 1106 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1107 1108 PetscFunctionBegin; 1109 if (op == MAT_ROW_ORIENTED) aij->roworiented = PETSC_TRUE; 1110 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = PETSC_FALSE; 1111 else if (op == MAT_ROWS_SORTED || 1112 op == MAT_ROWS_UNSORTED || 1113 op == MAT_COLUMNS_SORTED || 1114 op == MAT_COLUMNS_UNSORTED || 1115 op == MAT_SYMMETRIC || 1116 op == MAT_STRUCTURALLY_SYMMETRIC || 1117 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1118 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1119 op == MAT_NEW_NONZERO_LOCATION_ERR || 1120 op == MAT_NO_NEW_DIAGONALS || 1121 op == MAT_YES_NEW_DIAGONALS || 1122 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1123 op == MAT_USE_HASH_TABLE) { 1124 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1125 } else { 1126 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1127 } 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNC__ 1132 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries_SeqDense" 1133 int MatZeroEntries_SeqDense(Mat A) 1134 { 1135 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1136 int ierr; 1137 1138 PetscFunctionBegin; 1139 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(Scalar));CHKERRQ(ierr); 1140 PetscFunctionReturn(0); 1141 } 1142 1143 #undef __FUNC__ 1144 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize_SeqDense" 1145 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1146 { 1147 PetscFunctionBegin; 1148 *bs = 1; 1149 PetscFunctionReturn(0); 1150 } 1151 1152 #undef __FUNC__ 1153 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows_SeqDense" 1154 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1155 { 1156 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1157 int n = A->n,i,j,ierr,N,*rows; 1158 Scalar *slot; 1159 1160 PetscFunctionBegin; 1161 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1162 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1163 for (i=0; i<N; i++) { 1164 slot = l->v + rows[i]; 1165 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1166 } 1167 if (diag) { 1168 for (i=0; i<N; i++) { 1169 slot = l->v + (n+1)*rows[i]; 1170 *slot = *diag; 1171 } 1172 } 1173 ISRestoreIndices(is,&rows); 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNC__ 1178 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange_SeqDense" 1179 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1180 { 1181 PetscFunctionBegin; 1182 if (m) *m = 0; 1183 if (n) *n = A->m; 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNC__ 1188 #define __FUNC__ /*<a name=""></a>*/"MatGetArray_SeqDense" 1189 int MatGetArray_SeqDense(Mat A,Scalar **array) 1190 { 1191 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1192 1193 PetscFunctionBegin; 1194 *array = mat->v; 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNC__ 1199 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray_SeqDense" 1200 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1201 { 1202 PetscFunctionBegin; 1203 *array = 0; /* user cannot accidently use the array later */ 1204 PetscFunctionReturn(0); 1205 } 1206 1207 #undef __FUNC__ 1208 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix_SeqDense" 1209 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1210 { 1211 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1212 int i,j,ierr,m = A->m,*irow,*icol,nrows,ncols; 1213 Scalar *av,*bv,*v = mat->v; 1214 Mat newmat; 1215 1216 PetscFunctionBegin; 1217 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1218 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1219 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1220 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1221 1222 /* Check submatrixcall */ 1223 if (scall == MAT_REUSE_MATRIX) { 1224 int n_cols,n_rows; 1225 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1226 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1227 newmat = *B; 1228 } else { 1229 /* Create and fill new matrix */ 1230 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1231 } 1232 1233 /* Now extract the data pointers and do the copy,column at a time */ 1234 bv = ((Mat_SeqDense*)newmat->data)->v; 1235 1236 for (i=0; i<ncols; i++) { 1237 av = v + m*icol[i]; 1238 for (j=0; j<nrows; j++) { 1239 *bv++ = av[irow[j]]; 1240 } 1241 } 1242 1243 /* Assemble the matrices so that the correct flags are set */ 1244 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1245 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1246 1247 /* Free work space */ 1248 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1249 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1250 *B = newmat; 1251 PetscFunctionReturn(0); 1252 } 1253 1254 #undef __FUNC__ 1255 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices_SeqDense" 1256 int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1257 { 1258 int ierr,i; 1259 1260 PetscFunctionBegin; 1261 if (scall == MAT_INITIAL_MATRIX) { 1262 *B = (Mat*)PetscMalloc((n+1)*sizeof(Mat));CHKPTRQ(*B); 1263 } 1264 1265 for (i=0; i<n; i++) { 1266 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1267 } 1268 PetscFunctionReturn(0); 1269 } 1270 1271 #undef __FUNC__ 1272 #define __FUNC__ /*<a name=""></a>*/"MatCopy_SeqDense" 1273 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1274 { 1275 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1276 int ierr; 1277 PetscTruth flag; 1278 1279 PetscFunctionBegin; 1280 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flag);CHKERRQ(ierr); 1281 if (!flag) { 1282 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1283 PetscFunctionReturn(0); 1284 } 1285 if (A->m != B->m || A->n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1286 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(Scalar));CHKERRQ(ierr); 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNC__ 1291 #define __FUNC__ /*<a name="MatSetUpPreallocation_SeqDense"></a>*/"MatSetUpPreallocation_SeqDense" 1292 int MatSetUpPreallocation_SeqDense(Mat A) 1293 { 1294 int ierr; 1295 1296 PetscFunctionBegin; 1297 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /* -------------------------------------------------------------------*/ 1302 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1303 MatGetRow_SeqDense, 1304 MatRestoreRow_SeqDense, 1305 MatMult_SeqDense, 1306 MatMultAdd_SeqDense, 1307 MatMultTranspose_SeqDense, 1308 MatMultTransposeAdd_SeqDense, 1309 MatSolve_SeqDense, 1310 MatSolveAdd_SeqDense, 1311 MatSolveTranspose_SeqDense, 1312 MatSolveTransposeAdd_SeqDense, 1313 MatLUFactor_SeqDense, 1314 MatCholeskyFactor_SeqDense, 1315 MatRelax_SeqDense, 1316 MatTranspose_SeqDense, 1317 MatGetInfo_SeqDense, 1318 MatEqual_SeqDense, 1319 MatGetDiagonal_SeqDense, 1320 MatDiagonalScale_SeqDense, 1321 MatNorm_SeqDense, 1322 0, 1323 0, 1324 0, 1325 MatSetOption_SeqDense, 1326 MatZeroEntries_SeqDense, 1327 MatZeroRows_SeqDense, 1328 MatLUFactorSymbolic_SeqDense, 1329 MatLUFactorNumeric_SeqDense, 1330 MatCholeskyFactorSymbolic_SeqDense, 1331 MatCholeskyFactorNumeric_SeqDense, 1332 MatSetUpPreallocation_SeqDense, 1333 0, 1334 MatGetOwnershipRange_SeqDense, 1335 0, 1336 0, 1337 MatGetArray_SeqDense, 1338 MatRestoreArray_SeqDense, 1339 MatDuplicate_SeqDense, 1340 0, 1341 0, 1342 0, 1343 0, 1344 MatAXPY_SeqDense, 1345 MatGetSubMatrices_SeqDense, 1346 0, 1347 MatGetValues_SeqDense, 1348 MatCopy_SeqDense, 1349 0, 1350 MatScale_SeqDense, 1351 0, 1352 0, 1353 0, 1354 MatGetBlockSize_SeqDense, 1355 0, 1356 0, 1357 0, 1358 0, 1359 0, 1360 0, 1361 0, 1362 0, 1363 0, 1364 0, 1365 MatDestroy_SeqDense, 1366 MatView_SeqDense, 1367 MatGetMaps_Petsc}; 1368 1369 #undef __FUNC__ 1370 #define __FUNC__ /*<a name=""></a>*/"MatCreateSeqDense" 1371 /*@C 1372 MatCreateSeqDense - Creates a sequential dense matrix that 1373 is stored in column major order (the usual Fortran 77 manner). Many 1374 of the matrix operations use the BLAS and LAPACK routines. 1375 1376 Collective on MPI_Comm 1377 1378 Input Parameters: 1379 + comm - MPI communicator, set to PETSC_COMM_SELF 1380 . m - number of rows 1381 . n - number of columns 1382 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1383 to control all matrix memory allocation. 1384 1385 Output Parameter: 1386 . A - the matrix 1387 1388 Notes: 1389 The data input variable is intended primarily for Fortran programmers 1390 who wish to allocate their own matrix memory space. Most users should 1391 set data=PETSC_NULL. 1392 1393 Level: intermediate 1394 1395 .keywords: dense, matrix, LAPACK, BLAS 1396 1397 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1398 @*/ 1399 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1400 { 1401 int ierr; 1402 1403 PetscFunctionBegin; 1404 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1405 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1406 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1407 PetscFunctionReturn(0); 1408 } 1409 1410 #undef __FUNC__ 1411 #define __FUNC__ /*<a name=""></a>*/"MatSeqDensePreallocation" 1412 /*@C 1413 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1414 1415 Collective on MPI_Comm 1416 1417 Input Parameters: 1418 + A - the matrix 1419 - data - the array (or PETSC_NULL) 1420 1421 Notes: 1422 The data input variable is intended primarily for Fortran programmers 1423 who wish to allocate their own matrix memory space. Most users should 1424 set data=PETSC_NULL. 1425 1426 Level: intermediate 1427 1428 .keywords: dense, matrix, LAPACK, BLAS 1429 1430 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1431 @*/ 1432 int MatSeqDenseSetPreallocation(Mat B,Scalar *data) 1433 { 1434 Mat_SeqDense *b; 1435 int ierr; 1436 PetscTruth flg2; 1437 1438 PetscFunctionBegin; 1439 ierr = PetscTypeCompare((PetscObject)B,MATSEQDENSE,&flg2);CHKERRQ(ierr); 1440 if (!flg2) PetscFunctionReturn(0); 1441 B->preallocated = PETSC_TRUE; 1442 b = (Mat_SeqDense*)B->data; 1443 if (!data) { 1444 b->v = (Scalar*)PetscMalloc((B->m*B->n+1)*sizeof(Scalar));CHKPTRQ(b->v); 1445 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(Scalar));CHKERRQ(ierr); 1446 b->user_alloc = PETSC_FALSE; 1447 PLogObjectMemory(B,B->n*B->m*sizeof(Scalar)); 1448 } else { /* user-allocated storage */ 1449 b->v = data; 1450 b->user_alloc = PETSC_TRUE; 1451 } 1452 PetscFunctionReturn(0); 1453 } 1454 1455 EXTERN_C_BEGIN 1456 #undef __FUNC__ 1457 #define __FUNC__ /*<a name=""></a>*/"MatCreate_SeqDense" 1458 int MatCreate_SeqDense(Mat B) 1459 { 1460 Mat_SeqDense *b; 1461 int ierr,size; 1462 1463 PetscFunctionBegin; 1464 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1465 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1466 1467 B->m = B->M = PetscMax(B->m,B->M); 1468 B->n = B->N = PetscMax(B->n,B->N); 1469 1470 b = (Mat_SeqDense*)PetscMalloc(sizeof(Mat_SeqDense));CHKPTRQ(b); 1471 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1472 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1473 B->factor = 0; 1474 B->mapping = 0; 1475 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1476 B->data = (void*)b; 1477 1478 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1479 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1480 1481 b->pivots = 0; 1482 b->roworiented = PETSC_TRUE; 1483 b->v = 0; 1484 1485 PetscFunctionReturn(0); 1486 } 1487 EXTERN_C_END 1488 1489 1490 1491 1492