1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.167 1999/04/19 22:11:51 bsmith Exp balay $"; 3 #endif 4 /* 5 Defines the basic matrix operations for sequential dense. 6 */ 7 8 #include "src/mat/impls/dense/seq/dense.h" 9 #include "pinclude/blaslapack.h" 10 11 #undef __FUNC__ 12 #define __FUNC__ "MatAXPY_SeqDense" 13 int MatAXPY_SeqDense(Scalar *alpha,Mat X,Mat Y) 14 { 15 Mat_SeqDense *x = (Mat_SeqDense*) X->data,*y = (Mat_SeqDense*) Y->data; 16 int N = x->m*x->n, one = 1; 17 18 PetscFunctionBegin; 19 BLaxpy_( &N, alpha, x->v, &one, y->v, &one ); 20 PLogFlops(2*N-1); 21 PetscFunctionReturn(0); 22 } 23 24 #undef __FUNC__ 25 #define __FUNC__ "MatGetInfo_SeqDense" 26 int MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 27 { 28 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 29 int i,N = a->m*a->n,count = 0; 30 Scalar *v = a->v; 31 32 PetscFunctionBegin; 33 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 34 35 info->rows_global = (double)a->m; 36 info->columns_global = (double)a->n; 37 info->rows_local = (double)a->m; 38 info->columns_local = (double)a->n; 39 info->block_size = 1.0; 40 info->nz_allocated = (double)N; 41 info->nz_used = (double)count; 42 info->nz_unneeded = (double)(N-count); 43 info->assemblies = (double)A->num_ass; 44 info->mallocs = 0; 45 info->memory = A->mem; 46 info->fill_ratio_given = 0; 47 info->fill_ratio_needed = 0; 48 info->factor_mallocs = 0; 49 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNC__ 54 #define __FUNC__ "MatScale_SeqDense" 55 int MatScale_SeqDense(Scalar *alpha,Mat inA) 56 { 57 Mat_SeqDense *a = (Mat_SeqDense *) inA->data; 58 int one = 1, nz; 59 60 PetscFunctionBegin; 61 nz = a->m*a->n; 62 BLscal_( &nz, alpha, a->v, &one ); 63 PLogFlops(nz); 64 PetscFunctionReturn(0); 65 } 66 67 /* ---------------------------------------------------------------*/ 68 /* COMMENT: I have chosen to hide column permutation in the pivots, 69 rather than put it in the Mat->col slot.*/ 70 #undef __FUNC__ 71 #define __FUNC__ "MatLUFactor_SeqDense" 72 int MatLUFactor_SeqDense(Mat A,IS row,IS col,double f) 73 { 74 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 75 int info; 76 77 PetscFunctionBegin; 78 if (!mat->pivots) { 79 mat->pivots = (int *) PetscMalloc((mat->m+1)*sizeof(int));CHKPTRQ(mat->pivots); 80 PLogObjectMemory(A,mat->m*sizeof(int)); 81 } 82 if (!mat->m || !mat->n) PetscFunctionReturn(0); 83 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 84 if (info<0) SETERRQ(PETSC_ERR_LIB,0,"Bad argument to LU factorization"); 85 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Bad LU factorization"); 86 A->factor = FACTOR_LU; 87 PLogFlops((2*mat->n*mat->n*mat->n)/3); 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNC__ 92 #define __FUNC__ "MatDuplicate_SeqDense" 93 int MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 94 { 95 Mat_SeqDense *mat = (Mat_SeqDense *) A->data, *l; 96 int ierr; 97 Mat newi; 98 99 PetscFunctionBegin; 100 ierr = MatCreateSeqDense(A->comm,mat->m,mat->n,PETSC_NULL,&newi);CHKERRQ(ierr); 101 l = (Mat_SeqDense *) newi->data; 102 if (cpvalues == MAT_COPY_VALUES) { 103 ierr = PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar));CHKERRQ(ierr); 104 } 105 newi->assembled = PETSC_TRUE; 106 *newmat = newi; 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNC__ 111 #define __FUNC__ "MatLUFactorSymbolic_SeqDense" 112 int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,double f,Mat *fact) 113 { 114 int ierr; 115 116 PetscFunctionBegin; 117 ierr = MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr); 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNC__ 122 #define __FUNC__ "MatLUFactorNumeric_SeqDense" 123 int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 124 { 125 Mat_SeqDense *mat = (Mat_SeqDense*) A->data, *l = (Mat_SeqDense*) (*fact)->data; 126 int ierr; 127 128 PetscFunctionBegin; 129 /* copy the numerical values */ 130 ierr = PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar));CHKERRQ(ierr); 131 (*fact)->factor = 0; 132 ierr = MatLUFactor(*fact,0,0,1.0);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNC__ 137 #define __FUNC__ "MatCholeskyFactorSymbolic_SeqDense" 138 int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,double f,Mat *fact) 139 { 140 int ierr; 141 142 PetscFunctionBegin; 143 ierr = MatConvert(A,MATSAME,fact);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNC__ 148 #define __FUNC__ "MatCholeskyFactorNumeric_SeqDense" 149 int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 150 { 151 int ierr; 152 153 PetscFunctionBegin; 154 ierr = MatCholeskyFactor(*fact,0,1.0);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNC__ 159 #define __FUNC__ "MatCholeskyFactor_SeqDense" 160 int MatCholeskyFactor_SeqDense(Mat A,IS perm,double f) 161 { 162 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 163 int info; 164 165 PetscFunctionBegin; 166 if (mat->pivots) { 167 PetscFree(mat->pivots); 168 PLogObjectMemory(A,-mat->m*sizeof(int)); 169 mat->pivots = 0; 170 } 171 if (!mat->m || !mat->n) PetscFunctionReturn(0); 172 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 173 if (info) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,0,"Bad factorization"); 174 A->factor = FACTOR_CHOLESKY; 175 PLogFlops((mat->n*mat->n*mat->n)/3); 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNC__ 180 #define __FUNC__ "MatSolve_SeqDense" 181 int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 182 { 183 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 184 int one = 1, info, ierr; 185 Scalar *x, *y; 186 187 PetscFunctionBegin; 188 if (!mat->m || !mat->n) PetscFunctionReturn(0); 189 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 190 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 191 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 192 if (A->factor == FACTOR_LU) { 193 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 194 } else if (A->factor == FACTOR_CHOLESKY){ 195 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 196 } 197 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Matrix must be factored to solve"); 198 if (info) SETERRQ(PETSC_ERR_LIB,0,"MBad solve"); 199 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 200 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 201 PLogFlops(mat->n*mat->n - mat->n); 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNC__ 206 #define __FUNC__ "MatSolveTrans_SeqDense" 207 int MatSolveTrans_SeqDense(Mat A,Vec xx,Vec yy) 208 { 209 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 210 int ierr,one = 1, info; 211 Scalar *x, *y; 212 213 PetscFunctionBegin; 214 if (!mat->m || !mat->n) PetscFunctionReturn(0); 215 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 216 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 217 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 218 /* assume if pivots exist then use LU; else Cholesky */ 219 if (mat->pivots) { 220 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 221 } else { 222 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 223 } 224 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 225 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 226 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 227 PLogFlops(mat->n*mat->n - mat->n); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNC__ 232 #define __FUNC__ "MatSolveAdd_SeqDense" 233 int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 234 { 235 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 236 int ierr,one = 1, info; 237 Scalar *x, *y, sone = 1.0; 238 Vec tmp = 0; 239 240 PetscFunctionBegin; 241 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 242 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 243 if (!mat->m || !mat->n) PetscFunctionReturn(0); 244 if (yy == zz) { 245 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 246 PLogObjectParent(A,tmp); 247 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 248 } 249 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 250 /* assume if pivots exist then use LU; else Cholesky */ 251 if (mat->pivots) { 252 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 253 } else { 254 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 255 } 256 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 257 if (tmp) {ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 258 else {ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr);} 259 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 260 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 261 PLogFlops(mat->n*mat->n - mat->n); 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNC__ 266 #define __FUNC__ "MatSolveTransAdd_SeqDense" 267 int MatSolveTransAdd_SeqDense(Mat A,Vec xx,Vec zz, Vec yy) 268 { 269 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 270 int one = 1, info,ierr; 271 Scalar *x, *y, sone = 1.0; 272 Vec tmp; 273 274 PetscFunctionBegin; 275 if (!mat->m || !mat->n) PetscFunctionReturn(0); 276 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 277 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 278 if (yy == zz) { 279 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 280 PLogObjectParent(A,tmp); 281 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 282 } 283 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 284 /* assume if pivots exist then use LU; else Cholesky */ 285 if (mat->pivots) { 286 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 287 } else { 288 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 289 } 290 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 291 if (tmp) { 292 ierr = VecAXPY(&sone,tmp,yy);CHKERRQ(ierr); 293 ierr = VecDestroy(tmp);CHKERRQ(ierr); 294 } else { 295 ierr = VecAXPY(&sone,zz,yy);CHKERRQ(ierr); 296 } 297 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 298 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 299 PLogFlops(mat->n*mat->n - mat->n); 300 PetscFunctionReturn(0); 301 } 302 /* ------------------------------------------------------------------*/ 303 #undef __FUNC__ 304 #define __FUNC__ "MatRelax_SeqDense" 305 int MatRelax_SeqDense(Mat A,Vec bb,double omega,MatSORType flag, 306 double shift,int its,Vec xx) 307 { 308 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 309 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 310 int ierr, m = mat->m, i; 311 #if !defined(USE_PETSC_COMPLEX) 312 int o = 1; 313 #endif 314 315 PetscFunctionBegin; 316 if (flag & SOR_ZERO_INITIAL_GUESS) { 317 /* this is a hack fix, should have another version without the second BLdot */ 318 ierr = VecSet(&zero,xx);CHKERRQ(ierr); 319 } 320 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 321 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 322 while (its--) { 323 if (flag & SOR_FORWARD_SWEEP){ 324 for ( i=0; i<m; i++ ) { 325 #if defined(USE_PETSC_COMPLEX) 326 /* cannot use BLAS dot for complex because compiler/linker is 327 not happy about returning a double complex */ 328 int _i; 329 Scalar sum = b[i]; 330 for ( _i=0; _i<m; _i++ ) { 331 sum -= PetscConj(v[i+_i*m])*x[_i]; 332 } 333 xt = sum; 334 #else 335 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 336 #endif 337 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 338 } 339 } 340 if (flag & SOR_BACKWARD_SWEEP) { 341 for ( i=m-1; i>=0; i-- ) { 342 #if defined(USE_PETSC_COMPLEX) 343 /* cannot use BLAS dot for complex because compiler/linker is 344 not happy about returning a double complex */ 345 int _i; 346 Scalar sum = b[i]; 347 for ( _i=0; _i<m; _i++ ) { 348 sum -= PetscConj(v[i+_i*m])*x[_i]; 349 } 350 xt = sum; 351 #else 352 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 353 #endif 354 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 355 } 356 } 357 } 358 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 359 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 360 PetscFunctionReturn(0); 361 } 362 363 /* -----------------------------------------------------------------*/ 364 #undef __FUNC__ 365 #define __FUNC__ "MatMultTrans_SeqDense" 366 int MatMultTrans_SeqDense(Mat A,Vec xx,Vec yy) 367 { 368 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 369 Scalar *v = mat->v, *x, *y; 370 int ierr,_One=1;Scalar _DOne=1.0, _DZero=0.0; 371 372 PetscFunctionBegin; 373 if (!mat->m || !mat->n) PetscFunctionReturn(0); 374 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 375 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 376 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 377 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 378 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 379 PLogFlops(2*mat->m*mat->n - mat->n); 380 PetscFunctionReturn(0); 381 } 382 383 #undef __FUNC__ 384 #define __FUNC__ "MatMult_SeqDense" 385 int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 386 { 387 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 388 Scalar *v = mat->v, *x, *y; 389 int ierr,_One=1;Scalar _DOne=1.0, _DZero=0.0; 390 391 PetscFunctionBegin; 392 if (!mat->m || !mat->n) PetscFunctionReturn(0); 393 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 394 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 395 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DZero,y,&_One); 396 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 397 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 398 PLogFlops(2*mat->m*mat->n - mat->m); 399 PetscFunctionReturn(0); 400 } 401 402 #undef __FUNC__ 403 #define __FUNC__ "MatMultAdd_SeqDense" 404 int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 405 { 406 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 407 Scalar *v = mat->v, *x, *y; 408 int ierr,_One=1; Scalar _DOne=1.0; 409 410 PetscFunctionBegin; 411 if (!mat->m || !mat->n) PetscFunctionReturn(0); 412 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 413 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 414 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 415 LAgemv_( "N", &(mat->m), &(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 416 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 417 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 418 PLogFlops(2*mat->m*mat->n); 419 PetscFunctionReturn(0); 420 } 421 422 #undef __FUNC__ 423 #define __FUNC__ "MatMultTransAdd_SeqDense" 424 int MatMultTransAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 425 { 426 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 427 Scalar *v = mat->v, *x, *y; 428 int ierr,_One=1; 429 Scalar _DOne=1.0; 430 431 PetscFunctionBegin; 432 if (!mat->m || !mat->n) PetscFunctionReturn(0); 433 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 434 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 435 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 436 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DOne,y,&_One); 437 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 438 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 439 PLogFlops(2*mat->m*mat->n); 440 PetscFunctionReturn(0); 441 } 442 443 /* -----------------------------------------------------------------*/ 444 #undef __FUNC__ 445 #define __FUNC__ "MatGetRow_SeqDense" 446 int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 447 { 448 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 449 Scalar *v; 450 int i; 451 452 PetscFunctionBegin; 453 *ncols = mat->n; 454 if (cols) { 455 *cols = (int *) PetscMalloc((mat->n+1)*sizeof(int));CHKPTRQ(*cols); 456 for ( i=0; i<mat->n; i++ ) (*cols)[i] = i; 457 } 458 if (vals) { 459 *vals = (Scalar *) PetscMalloc((mat->n+1)*sizeof(Scalar));CHKPTRQ(*vals); 460 v = mat->v + row; 461 for ( i=0; i<mat->n; i++ ) {(*vals)[i] = *v; v += mat->m;} 462 } 463 PetscFunctionReturn(0); 464 } 465 466 #undef __FUNC__ 467 #define __FUNC__ "MatRestoreRow_SeqDense" 468 int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 469 { 470 if (cols) { PetscFree(*cols); } 471 if (vals) { PetscFree(*vals); } 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(USE_PETSC_BOPT_g) 489 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 490 #endif 491 for ( i=0; i<m; i++ ) { 492 if (indexm[i] < 0) {v++; continue;} 493 #if defined(USE_PETSC_BOPT_g) 494 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 495 #endif 496 mat->v[indexn[j]*mat->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(USE_PETSC_BOPT_g) 503 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 504 #endif 505 for ( i=0; i<m; i++ ) { 506 if (indexm[i] < 0) {v++; continue;} 507 #if defined(USE_PETSC_BOPT_g) 508 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 509 #endif 510 mat->v[indexn[j]*mat->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(USE_PETSC_BOPT_g) 519 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 520 #endif 521 for ( j=0; j<n; j++ ) { 522 if (indexn[j] < 0) { v++; continue;} 523 #if defined(USE_PETSC_BOPT_g) 524 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 525 #endif 526 mat->v[indexn[j]*mat->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(USE_PETSC_BOPT_g) 533 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 534 #endif 535 for ( j=0; j<n; j++ ) { 536 if (indexn[j] < 0) { v++; continue;} 537 #if defined(USE_PETSC_BOPT_g) 538 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 539 #endif 540 mat->v[indexn[j]*mat->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]*mat->m + indexm[i]]; 561 } 562 } 563 PetscFunctionReturn(0); 564 } 565 566 /* -----------------------------------------------------------------*/ 567 568 #include "sys.h" 569 570 #undef __FUNC__ 571 #define __FUNC__ "MatLoad_SeqDense" 572 int MatLoad_SeqDense(Viewer viewer,MatType type,Mat *A) 573 { 574 Mat_SeqDense *a; 575 Mat B; 576 int *scols, i, j, nz, ierr, fd, header[4], size; 577 int *rowlengths = 0, M, N, *cols; 578 Scalar *vals, *svals, *v,*w; 579 MPI_Comm comm = ((PetscObject)viewer)->comm; 580 581 PetscFunctionBegin; 582 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 583 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 584 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 585 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 586 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Not matrix object"); 587 M = header[1]; N = header[2]; nz = header[3]; 588 589 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 590 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 591 B = *A; 592 a = (Mat_SeqDense *) B->data; 593 v = a->v; 594 /* Allocate some temp space to read in the values and then flip them 595 from row major to column major */ 596 w = (Scalar *) PetscMalloc((M*N+1)*sizeof(Scalar));CHKPTRQ(w); 597 /* read in nonzero values */ 598 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 599 /* now flip the values and store them in the matrix*/ 600 for ( j=0; j<N; j++ ) { 601 for ( i=0; i<M; i++ ) { 602 *v++ =w[i*N+j]; 603 } 604 } 605 PetscFree(w); 606 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 607 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 608 } else { 609 /* read row lengths */ 610 rowlengths = (int*) PetscMalloc( (M+1)*sizeof(int) );CHKPTRQ(rowlengths); 611 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 612 613 /* create our matrix */ 614 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 615 B = *A; 616 a = (Mat_SeqDense *) B->data; 617 v = a->v; 618 619 /* read column indices and nonzeros */ 620 cols = scols = (int *) PetscMalloc( (nz+1)*sizeof(int) );CHKPTRQ(cols); 621 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 622 vals = svals = (Scalar *) PetscMalloc( (nz+1)*sizeof(Scalar) );CHKPTRQ(vals); 623 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 624 625 /* insert into matrix */ 626 for ( i=0; i<M; i++ ) { 627 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 628 svals += rowlengths[i]; scols += rowlengths[i]; 629 } 630 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 631 632 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 633 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 634 } 635 PetscFunctionReturn(0); 636 } 637 638 #include "sys.h" 639 640 #undef __FUNC__ 641 #define __FUNC__ "MatView_SeqDense_ASCII" 642 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 643 { 644 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 645 int ierr, i, j, format; 646 FILE *fd; 647 char *outputname; 648 Scalar *v; 649 650 PetscFunctionBegin; 651 ierr = ViewerASCIIGetPointer(viewer,&fd);CHKERRQ(ierr); 652 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 653 ierr = ViewerGetFormat(viewer,&format); 654 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 655 PetscFunctionReturn(0); /* do nothing for now */ 656 } 657 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 658 for ( i=0; i<a->m; i++ ) { 659 v = a->v + i; 660 fprintf(fd,"row %d:",i); 661 for ( j=0; j<a->n; j++ ) { 662 #if defined(USE_PETSC_COMPLEX) 663 if (PetscReal(*v) != 0.0 && PetscImaginary(*v) != 0.0) fprintf(fd," %d %g + %g i",j,PetscReal(*v),PetscImaginary(*v)); 664 else if (PetscReal(*v)) fprintf(fd," %d %g ",j,PetscReal(*v)); 665 #else 666 if (*v) fprintf(fd," %d %g ",j,*v); 667 #endif 668 v += a->m; 669 } 670 fprintf(fd,"\n"); 671 } 672 } else { 673 #if defined(USE_PETSC_COMPLEX) 674 int allreal = 1; 675 /* determine if matrix has all real values */ 676 v = a->v; 677 for ( i=0; i<a->m*a->n; i++ ) { 678 if (PetscImaginary(v[i])) { allreal = 0; break ;} 679 } 680 #endif 681 for ( i=0; i<a->m; i++ ) { 682 v = a->v + i; 683 for ( j=0; j<a->n; j++ ) { 684 #if defined(USE_PETSC_COMPLEX) 685 if (allreal) { 686 fprintf(fd,"%6.4e ",PetscReal(*v)); v += a->m; 687 } else { 688 fprintf(fd,"%6.4e + %6.4e i ",PetscReal(*v),PetscImaginary(*v)); v += a->m; 689 } 690 #else 691 fprintf(fd,"%6.4e ",*v); v += a->m; 692 #endif 693 } 694 fprintf(fd,"\n"); 695 } 696 } 697 fflush(fd); 698 PetscFunctionReturn(0); 699 } 700 701 #undef __FUNC__ 702 #define __FUNC__ "MatView_SeqDense_Binary" 703 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 704 { 705 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 706 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 707 int format; 708 Scalar *v, *anonz,*vals; 709 710 PetscFunctionBegin; 711 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 712 713 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 714 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 715 /* store the matrix as a dense matrix */ 716 col_lens = (int *) PetscMalloc( 4*sizeof(int) );CHKPTRQ(col_lens); 717 col_lens[0] = MAT_COOKIE; 718 col_lens[1] = m; 719 col_lens[2] = n; 720 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 721 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 722 PetscFree(col_lens); 723 724 /* write out matrix, by rows */ 725 vals = (Scalar *) PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(vals); 726 v = a->v; 727 for ( i=0; i<m; i++ ) { 728 for ( j=0; j<n; j++ ) { 729 vals[i + j*m] = *v++; 730 } 731 } 732 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 733 PetscFree(vals); 734 } else { 735 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) );CHKPTRQ(col_lens); 736 col_lens[0] = MAT_COOKIE; 737 col_lens[1] = m; 738 col_lens[2] = n; 739 col_lens[3] = nz; 740 741 /* store lengths of each row and write (including header) to file */ 742 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 743 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 744 745 /* Possibly should write in smaller increments, not whole matrix at once? */ 746 /* store column indices (zero start index) */ 747 ict = 0; 748 for ( i=0; i<m; i++ ) { 749 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 750 } 751 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 752 PetscFree(col_lens); 753 754 /* store nonzero values */ 755 anonz = (Scalar *) PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(anonz); 756 ict = 0; 757 for ( i=0; i<m; i++ ) { 758 v = a->v + i; 759 for ( j=0; j<n; j++ ) { 760 anonz[ict++] = *v; v += a->m; 761 } 762 } 763 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 764 PetscFree(anonz); 765 } 766 PetscFunctionReturn(0); 767 } 768 769 #undef __FUNC__ 770 #define __FUNC__ "MatView_SeqDense" 771 int MatView_SeqDense(Mat A,Viewer viewer) 772 { 773 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 774 ViewerType vtype; 775 int ierr; 776 777 PetscFunctionBegin; 778 ierr = ViewerGetType(viewer,&vtype);CHKERRQ(ierr); 779 780 if (PetscTypeCompare(vtype,SOCKET_VIEWER)) { 781 ierr = ViewerSocketPutScalar_Private(viewer,a->m,a->n,a->v);CHKERRQ(ierr); 782 } else if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 783 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 784 } else if (PetscTypeCompare(vtype,BINARY_VIEWER)) { 785 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 786 } else { 787 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 788 } 789 PetscFunctionReturn(0); 790 } 791 792 #undef __FUNC__ 793 #define __FUNC__ "MatDestroy_SeqDense" 794 int MatDestroy_SeqDense(Mat mat) 795 { 796 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 797 int ierr; 798 799 PetscFunctionBegin; 800 if (--mat->refct > 0) PetscFunctionReturn(0); 801 802 if (mat->mapping) { 803 ierr = ISLocalToGlobalMappingDestroy(mat->mapping);CHKERRQ(ierr); 804 } 805 if (mat->bmapping) { 806 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping);CHKERRQ(ierr); 807 } 808 #if defined(USE_PETSC_LOG) 809 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 810 #endif 811 if (l->pivots) PetscFree(l->pivots); 812 if (!l->user_alloc) PetscFree(l->v); 813 PetscFree(l); 814 if (mat->rmap) { 815 ierr = MapDestroy(mat->rmap);CHKERRQ(ierr); 816 } 817 if (mat->cmap) { 818 ierr = MapDestroy(mat->cmap);CHKERRQ(ierr); 819 } 820 PLogObjectDestroy(mat); 821 PetscHeaderDestroy(mat); 822 PetscFunctionReturn(0); 823 } 824 825 #undef __FUNC__ 826 #define __FUNC__ "MatTranspose_SeqDense" 827 int MatTranspose_SeqDense(Mat A,Mat *matout) 828 { 829 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 830 int k, j, m, n, ierr; 831 Scalar *v, tmp; 832 833 PetscFunctionBegin; 834 v = mat->v; m = mat->m; n = mat->n; 835 if (matout == PETSC_NULL) { /* in place transpose */ 836 if (m != n) { /* malloc temp to hold transpose */ 837 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 838 for ( j=0; j<m; j++ ) { 839 for ( k=0; k<n; k++ ) { 840 w[k + j*n] = v[j + k*m]; 841 } 842 } 843 ierr = PetscMemcpy(v,w,m*n*sizeof(Scalar));CHKERRQ(ierr); 844 PetscFree(w); 845 } else { 846 for ( j=0; j<m; j++ ) { 847 for ( k=0; k<j; k++ ) { 848 tmp = v[j + k*n]; 849 v[j + k*n] = v[k + j*n]; 850 v[k + j*n] = tmp; 851 } 852 } 853 } 854 } else { /* out-of-place transpose */ 855 int ierr; 856 Mat tmat; 857 Mat_SeqDense *tmatd; 858 Scalar *v2; 859 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 860 tmatd = (Mat_SeqDense *) tmat->data; 861 v = mat->v; v2 = tmatd->v; 862 for ( j=0; j<n; j++ ) { 863 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 864 } 865 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 866 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 867 *matout = tmat; 868 } 869 PetscFunctionReturn(0); 870 } 871 872 #undef __FUNC__ 873 #define __FUNC__ "MatEqual_SeqDense" 874 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 875 { 876 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 877 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 878 int i; 879 Scalar *v1 = mat1->v, *v2 = mat2->v; 880 881 PetscFunctionBegin; 882 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 883 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 884 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 885 for ( i=0; i<mat1->m*mat1->n; i++ ) { 886 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 887 v1++; v2++; 888 } 889 *flg = PETSC_TRUE; 890 PetscFunctionReturn(0); 891 } 892 893 #undef __FUNC__ 894 #define __FUNC__ "MatGetDiagonal_SeqDense" 895 int MatGetDiagonal_SeqDense(Mat A,Vec v) 896 { 897 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 898 int ierr,i, n, len; 899 Scalar *x, zero = 0.0; 900 901 PetscFunctionBegin; 902 ierr = VecSet(&zero,v);CHKERRQ(ierr); 903 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 904 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 905 len = PetscMin(mat->m,mat->n); 906 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 907 for ( i=0; i<len; i++ ) { 908 x[i] = mat->v[i*mat->m + i]; 909 } 910 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 911 PetscFunctionReturn(0); 912 } 913 914 #undef __FUNC__ 915 #define __FUNC__ "MatDiagonalScale_SeqDense" 916 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 917 { 918 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 919 Scalar *l,*r,x,*v; 920 int ierr,i,j,m = mat->m, n = mat->n; 921 922 PetscFunctionBegin; 923 if (ll) { 924 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 925 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 926 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 927 for ( i=0; i<m; i++ ) { 928 x = l[i]; 929 v = mat->v + i; 930 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 931 } 932 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 933 PLogFlops(n*m); 934 } 935 if (rr) { 936 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 937 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 938 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 939 for ( i=0; i<n; i++ ) { 940 x = r[i]; 941 v = mat->v + i*m; 942 for ( j=0; j<m; j++ ) { (*v++) *= x;} 943 } 944 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 945 PLogFlops(n*m); 946 } 947 PetscFunctionReturn(0); 948 } 949 950 #undef __FUNC__ 951 #define __FUNC__ "MatNorm_SeqDense" 952 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 953 { 954 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 955 Scalar *v = mat->v; 956 double sum = 0.0; 957 int i, j; 958 959 PetscFunctionBegin; 960 if (type == NORM_FROBENIUS) { 961 for (i=0; i<mat->n*mat->m; i++ ) { 962 #if defined(USE_PETSC_COMPLEX) 963 sum += PetscReal(PetscConj(*v)*(*v)); v++; 964 #else 965 sum += (*v)*(*v); v++; 966 #endif 967 } 968 *norm = sqrt(sum); 969 PLogFlops(2*mat->n*mat->m); 970 } else if (type == NORM_1) { 971 *norm = 0.0; 972 for ( j=0; j<mat->n; j++ ) { 973 sum = 0.0; 974 for ( i=0; i<mat->m; i++ ) { 975 sum += PetscAbsScalar(*v); v++; 976 } 977 if (sum > *norm) *norm = sum; 978 } 979 PLogFlops(mat->n*mat->m); 980 } else if (type == NORM_INFINITY) { 981 *norm = 0.0; 982 for ( j=0; j<mat->m; j++ ) { 983 v = mat->v + j; 984 sum = 0.0; 985 for ( i=0; i<mat->n; i++ ) { 986 sum += PetscAbsScalar(*v); v += mat->m; 987 } 988 if (sum > *norm) *norm = sum; 989 } 990 PLogFlops(mat->n*mat->m); 991 } else { 992 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 993 } 994 PetscFunctionReturn(0); 995 } 996 997 #undef __FUNC__ 998 #define __FUNC__ "MatSetOption_SeqDense" 999 int MatSetOption_SeqDense(Mat A,MatOption op) 1000 { 1001 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 1002 1003 PetscFunctionBegin; 1004 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 1005 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1006 else if (op == MAT_ROWS_SORTED || 1007 op == MAT_ROWS_UNSORTED || 1008 op == MAT_COLUMNS_SORTED || 1009 op == MAT_COLUMNS_UNSORTED || 1010 op == MAT_SYMMETRIC || 1011 op == MAT_STRUCTURALLY_SYMMETRIC || 1012 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1013 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1014 op == MAT_NEW_NONZERO_LOCATION_ERR || 1015 op == MAT_NO_NEW_DIAGONALS || 1016 op == MAT_YES_NEW_DIAGONALS || 1017 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1018 op == MAT_USE_HASH_TABLE) 1019 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1020 else { 1021 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1022 } 1023 PetscFunctionReturn(0); 1024 } 1025 1026 #undef __FUNC__ 1027 #define __FUNC__ "MatZeroEntries_SeqDense" 1028 int MatZeroEntries_SeqDense(Mat A) 1029 { 1030 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1031 int ierr; 1032 1033 PetscFunctionBegin; 1034 ierr = PetscMemzero(l->v,l->m*l->n*sizeof(Scalar));CHKERRQ(ierr); 1035 PetscFunctionReturn(0); 1036 } 1037 1038 #undef __FUNC__ 1039 #define __FUNC__ "MatGetBlockSize_SeqDense" 1040 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1041 { 1042 PetscFunctionBegin; 1043 *bs = 1; 1044 PetscFunctionReturn(0); 1045 } 1046 1047 #undef __FUNC__ 1048 #define __FUNC__ "MatZeroRows_SeqDense" 1049 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1050 { 1051 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1052 int n = l->n, i, j,ierr,N, *rows; 1053 Scalar *slot; 1054 1055 PetscFunctionBegin; 1056 ierr = ISGetSize(is,&N);CHKERRQ(ierr); 1057 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1058 for ( i=0; i<N; i++ ) { 1059 slot = l->v + rows[i]; 1060 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1061 } 1062 if (diag) { 1063 for ( i=0; i<N; i++ ) { 1064 slot = l->v + (n+1)*rows[i]; 1065 *slot = *diag; 1066 } 1067 } 1068 ISRestoreIndices(is,&rows); 1069 PetscFunctionReturn(0); 1070 } 1071 1072 #undef __FUNC__ 1073 #define __FUNC__ "MatGetSize_SeqDense" 1074 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1075 { 1076 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1077 1078 PetscFunctionBegin; 1079 *m = mat->m; *n = mat->n; 1080 PetscFunctionReturn(0); 1081 } 1082 1083 #undef __FUNC__ 1084 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1085 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1086 { 1087 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1088 1089 PetscFunctionBegin; 1090 *m = 0; *n = mat->m; 1091 PetscFunctionReturn(0); 1092 } 1093 1094 #undef __FUNC__ 1095 #define __FUNC__ "MatGetArray_SeqDense" 1096 int MatGetArray_SeqDense(Mat A,Scalar **array) 1097 { 1098 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1099 1100 PetscFunctionBegin; 1101 *array = mat->v; 1102 PetscFunctionReturn(0); 1103 } 1104 1105 #undef __FUNC__ 1106 #define __FUNC__ "MatRestoreArray_SeqDense" 1107 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1108 { 1109 PetscFunctionBegin; 1110 *array = 0; /* user cannot accidently use the array later */ 1111 PetscFunctionReturn(0); 1112 } 1113 1114 #undef __FUNC__ 1115 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1116 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1117 { 1118 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1119 int i, j, ierr, m = mat->m, *irow, *icol, nrows, ncols; 1120 Scalar *av, *bv, *v = mat->v; 1121 Mat newmat; 1122 1123 PetscFunctionBegin; 1124 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1125 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1126 ierr = ISGetSize(isrow,&nrows);CHKERRQ(ierr); 1127 ierr = ISGetSize(iscol,&ncols);CHKERRQ(ierr); 1128 1129 /* Check submatrixcall */ 1130 if (scall == MAT_REUSE_MATRIX) { 1131 int n_cols,n_rows; 1132 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1133 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1134 newmat = *B; 1135 } else { 1136 /* Create and fill new matrix */ 1137 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1138 } 1139 1140 /* Now extract the data pointers and do the copy, column at a time */ 1141 bv = ((Mat_SeqDense*)newmat->data)->v; 1142 1143 for ( i=0; i<ncols; i++ ) { 1144 av = v + m*icol[i]; 1145 for (j=0; j<nrows; j++ ) { 1146 *bv++ = av[irow[j]]; 1147 } 1148 } 1149 1150 /* Assemble the matrices so that the correct flags are set */ 1151 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1152 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1153 1154 /* Free work space */ 1155 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1156 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1157 *B = newmat; 1158 PetscFunctionReturn(0); 1159 } 1160 1161 #undef __FUNC__ 1162 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1163 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B) 1164 { 1165 int ierr,i; 1166 1167 PetscFunctionBegin; 1168 if (scall == MAT_INITIAL_MATRIX) { 1169 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) );CHKPTRQ(*B); 1170 } 1171 1172 for ( i=0; i<n; i++ ) { 1173 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1174 } 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNC__ 1179 #define __FUNC__ "MatCopy_SeqDense" 1180 int MatCopy_SeqDense(Mat A, Mat B,MatStructure str) 1181 { 1182 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1183 int ierr; 1184 1185 PetscFunctionBegin; 1186 if (B->type != MATSEQDENSE) { 1187 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1188 PetscFunctionReturn(0); 1189 } 1190 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1191 ierr = PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar));CHKERRQ(ierr); 1192 PetscFunctionReturn(0); 1193 } 1194 1195 /* -------------------------------------------------------------------*/ 1196 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1197 MatGetRow_SeqDense, 1198 MatRestoreRow_SeqDense, 1199 MatMult_SeqDense, 1200 MatMultAdd_SeqDense, 1201 MatMultTrans_SeqDense, 1202 MatMultTransAdd_SeqDense, 1203 MatSolve_SeqDense, 1204 MatSolveAdd_SeqDense, 1205 MatSolveTrans_SeqDense, 1206 MatSolveTransAdd_SeqDense, 1207 MatLUFactor_SeqDense, 1208 MatCholeskyFactor_SeqDense, 1209 MatRelax_SeqDense, 1210 MatTranspose_SeqDense, 1211 MatGetInfo_SeqDense, 1212 MatEqual_SeqDense, 1213 MatGetDiagonal_SeqDense, 1214 MatDiagonalScale_SeqDense, 1215 MatNorm_SeqDense, 1216 0, 1217 0, 1218 0, 1219 MatSetOption_SeqDense, 1220 MatZeroEntries_SeqDense, 1221 MatZeroRows_SeqDense, 1222 MatLUFactorSymbolic_SeqDense, 1223 MatLUFactorNumeric_SeqDense, 1224 MatCholeskyFactorSymbolic_SeqDense, 1225 MatCholeskyFactorNumeric_SeqDense, 1226 MatGetSize_SeqDense, 1227 MatGetSize_SeqDense, 1228 MatGetOwnershipRange_SeqDense, 1229 0, 1230 0, 1231 MatGetArray_SeqDense, 1232 MatRestoreArray_SeqDense, 1233 MatDuplicate_SeqDense, 1234 0, 1235 0, 1236 0, 1237 0, 1238 MatAXPY_SeqDense, 1239 MatGetSubMatrices_SeqDense, 1240 0, 1241 MatGetValues_SeqDense, 1242 MatCopy_SeqDense, 1243 0, 1244 MatScale_SeqDense, 1245 0, 1246 0, 1247 0, 1248 MatGetBlockSize_SeqDense, 1249 0, 1250 0, 1251 0, 1252 0, 1253 0, 1254 0, 1255 0, 1256 0, 1257 0, 1258 0, 1259 0, 1260 0, 1261 MatGetMaps_Petsc}; 1262 1263 #undef __FUNC__ 1264 #define __FUNC__ "MatCreateSeqDense" 1265 /*@C 1266 MatCreateSeqDense - Creates a sequential dense matrix that 1267 is stored in column major order (the usual Fortran 77 manner). Many 1268 of the matrix operations use the BLAS and LAPACK routines. 1269 1270 Collective on MPI_Comm 1271 1272 Input Parameters: 1273 + comm - MPI communicator, set to PETSC_COMM_SELF 1274 . m - number of rows 1275 . n - number of columns 1276 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1277 to control all matrix memory allocation. 1278 1279 Output Parameter: 1280 . A - the matrix 1281 1282 Notes: 1283 The data input variable is intended primarily for Fortran programmers 1284 who wish to allocate their own matrix memory space. Most users should 1285 set data=PETSC_NULL. 1286 1287 Level: intermediate 1288 1289 .keywords: dense, matrix, LAPACK, BLAS 1290 1291 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1292 @*/ 1293 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1294 { 1295 Mat B; 1296 Mat_SeqDense *b; 1297 int ierr,flg,size; 1298 1299 PetscFunctionBegin; 1300 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1301 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1302 1303 *A = 0; 1304 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,"Mat",comm,MatDestroy,MatView); 1305 PLogObjectCreate(B); 1306 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense));CHKPTRQ(b); 1307 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1308 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1309 B->ops->destroy = MatDestroy_SeqDense; 1310 B->ops->view = MatView_SeqDense; 1311 B->factor = 0; 1312 B->mapping = 0; 1313 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1314 B->data = (void *) b; 1315 1316 b->m = m; B->m = m; B->M = m; 1317 b->n = n; B->n = n; B->N = n; 1318 1319 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1320 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1321 1322 b->pivots = 0; 1323 b->roworiented = 1; 1324 if (data == PETSC_NULL) { 1325 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(b->v); 1326 ierr = PetscMemzero(b->v,m*n*sizeof(Scalar));CHKERRQ(ierr); 1327 b->user_alloc = 0; 1328 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1329 } else { /* user-allocated storage */ 1330 b->v = data; 1331 b->user_alloc = 1; 1332 } 1333 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 1334 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr);} 1335 1336 *A = B; 1337 PetscFunctionReturn(0); 1338 } 1339 1340 1341 1342 1343 1344 1345