1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.158 1998/10/01 21:32:47 balay Exp bsmith $"; 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 #include "pinclude/pviewer.h" 11 12 #undef __FUNC__ 13 #define __FUNC__ "MatAXPY_SeqDense" 14 int MatAXPY_SeqDense(Scalar *alpha,Mat X,Mat Y) 15 { 16 Mat_SeqDense *x = (Mat_SeqDense*) X->data,*y = (Mat_SeqDense*) Y->data; 17 int N = x->m*x->n, one = 1; 18 19 PetscFunctionBegin; 20 BLaxpy_( &N, alpha, x->v, &one, y->v, &one ); 21 PLogFlops(2*N-1); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNC__ 26 #define __FUNC__ "MatGetInfo_SeqDense" 27 int MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 28 { 29 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 30 int i,N = a->m*a->n,count = 0; 31 Scalar *v = a->v; 32 33 PetscFunctionBegin; 34 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 35 36 info->rows_global = (double)a->m; 37 info->columns_global = (double)a->n; 38 info->rows_local = (double)a->m; 39 info->columns_local = (double)a->n; 40 info->block_size = 1.0; 41 info->nz_allocated = (double)N; 42 info->nz_used = (double)count; 43 info->nz_unneeded = (double)(N-count); 44 info->assemblies = (double)A->num_ass; 45 info->mallocs = 0; 46 info->memory = A->mem; 47 info->fill_ratio_given = 0; 48 info->fill_ratio_needed = 0; 49 info->factor_mallocs = 0; 50 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNC__ 55 #define __FUNC__ "MatScale_SeqDense" 56 int MatScale_SeqDense(Scalar *alpha,Mat inA) 57 { 58 Mat_SeqDense *a = (Mat_SeqDense *) inA->data; 59 int one = 1, nz; 60 61 PetscFunctionBegin; 62 nz = a->m*a->n; 63 BLscal_( &nz, alpha, a->v, &one ); 64 PLogFlops(nz); 65 PetscFunctionReturn(0); 66 } 67 68 /* ---------------------------------------------------------------*/ 69 /* COMMENT: I have chosen to hide column permutation in the pivots, 70 rather than put it in the Mat->col slot.*/ 71 #undef __FUNC__ 72 #define __FUNC__ "MatLUFactor_SeqDense" 73 int MatLUFactor_SeqDense(Mat A,IS row,IS col,double f) 74 { 75 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 76 int info; 77 78 PetscFunctionBegin; 79 if (!mat->pivots) { 80 mat->pivots = (int *) PetscMalloc((mat->m+1)*sizeof(int));CHKPTRQ(mat->pivots); 81 PLogObjectMemory(A,mat->m*sizeof(int)); 82 } 83 if (!mat->m || !mat->n) PetscFunctionReturn(0); 84 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 85 if (info<0) SETERRQ(PETSC_ERR_LIB,0,"Bad argument to LU factorization"); 86 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Bad LU factorization"); 87 A->factor = FACTOR_LU; 88 PLogFlops((2*mat->n*mat->n*mat->n)/3); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNC__ 93 #define __FUNC__ "MatDuplicate_SeqDense" 94 int MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 95 { 96 Mat_SeqDense *mat = (Mat_SeqDense *) A->data, *l; 97 int ierr; 98 Mat newi; 99 100 PetscFunctionBegin; 101 ierr = MatCreateSeqDense(A->comm,mat->m,mat->n,PETSC_NULL,&newi); CHKERRQ(ierr); 102 l = (Mat_SeqDense *) newi->data; 103 if (cpvalues == MAT_COPY_VALUES) { 104 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 105 } 106 newi->assembled = PETSC_TRUE; 107 *newmat = newi; 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNC__ 112 #define __FUNC__ "MatLUFactorSymbolic_SeqDense" 113 int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,double f,Mat *fact) 114 { 115 int ierr; 116 117 PetscFunctionBegin; 118 ierr = MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);CHKERRQ(ierr); 119 PetscFunctionReturn(0); 120 } 121 122 #undef __FUNC__ 123 #define __FUNC__ "MatLUFactorNumeric_SeqDense" 124 int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 125 { 126 Mat_SeqDense *mat = (Mat_SeqDense*) A->data, *l = (Mat_SeqDense*) (*fact)->data; 127 int ierr; 128 129 PetscFunctionBegin; 130 /* copy the numerical values */ 131 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 132 (*fact)->factor = 0; 133 ierr = MatLUFactor(*fact,0,0,1.0);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNC__ 138 #define __FUNC__ "MatCholeskyFactorSymbolic_SeqDense" 139 int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,double f,Mat *fact) 140 { 141 int ierr; 142 143 PetscFunctionBegin; 144 ierr = MatConvert(A,MATSAME,fact);CHKERRQ(ierr); 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNC__ 149 #define __FUNC__ "MatCholeskyFactorNumeric_SeqDense" 150 int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 151 { 152 int ierr; 153 154 PetscFunctionBegin; 155 ierr = MatCholeskyFactor(*fact,0,1.0);CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 159 #undef __FUNC__ 160 #define __FUNC__ "MatCholeskyFactor_SeqDense" 161 int MatCholeskyFactor_SeqDense(Mat A,IS perm,double f) 162 { 163 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 164 int info; 165 166 PetscFunctionBegin; 167 if (mat->pivots) { 168 PetscFree(mat->pivots); 169 PLogObjectMemory(A,-mat->m*sizeof(int)); 170 mat->pivots = 0; 171 } 172 if (!mat->m || !mat->n) PetscFunctionReturn(0); 173 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 174 if (info) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,0,"Bad factorization"); 175 A->factor = FACTOR_CHOLESKY; 176 PLogFlops((mat->n*mat->n*mat->n)/3); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNC__ 181 #define __FUNC__ "MatSolve_SeqDense" 182 int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 183 { 184 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 185 int one = 1, info, ierr; 186 Scalar *x, *y; 187 188 PetscFunctionBegin; 189 if (!mat->m || !mat->n) PetscFunctionReturn(0); 190 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 191 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 192 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 193 if (A->factor == FACTOR_LU) { 194 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 195 } else if (A->factor == FACTOR_CHOLESKY){ 196 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 197 } 198 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Matrix must be factored to solve"); 199 if (info) SETERRQ(PETSC_ERR_LIB,0,"MBad solve"); 200 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 201 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 202 PLogFlops(mat->n*mat->n - mat->n); 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNC__ 207 #define __FUNC__ "MatSolveTrans_SeqDense" 208 int MatSolveTrans_SeqDense(Mat A,Vec xx,Vec yy) 209 { 210 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 211 int ierr,one = 1, info; 212 Scalar *x, *y; 213 214 PetscFunctionBegin; 215 if (!mat->m || !mat->n) PetscFunctionReturn(0); 216 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 217 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 218 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 219 /* assume if pivots exist then use LU; else Cholesky */ 220 if (mat->pivots) { 221 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 222 } else { 223 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 224 } 225 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 226 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 227 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 228 PLogFlops(mat->n*mat->n - mat->n); 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNC__ 233 #define __FUNC__ "MatSolveAdd_SeqDense" 234 int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 235 { 236 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 237 int ierr,one = 1, info; 238 Scalar *x, *y, sone = 1.0; 239 Vec tmp = 0; 240 241 PetscFunctionBegin; 242 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 243 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 244 if (!mat->m || !mat->n) PetscFunctionReturn(0); 245 if (yy == zz) { 246 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 247 PLogObjectParent(A,tmp); 248 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 249 } 250 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 251 /* assume if pivots exist then use LU; else Cholesky */ 252 if (mat->pivots) { 253 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 254 } else { 255 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 256 } 257 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 258 if (tmp) {ierr = VecAXPY(&sone,tmp,yy); CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 259 else {ierr = VecAXPY(&sone,zz,yy); CHKERRQ(ierr);} 260 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 261 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 262 PLogFlops(mat->n*mat->n - mat->n); 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNC__ 267 #define __FUNC__ "MatSolveTransAdd_SeqDense" 268 int MatSolveTransAdd_SeqDense(Mat A,Vec xx,Vec zz, Vec yy) 269 { 270 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 271 int one = 1, info,ierr; 272 Scalar *x, *y, sone = 1.0; 273 Vec tmp; 274 275 PetscFunctionBegin; 276 if (!mat->m || !mat->n) PetscFunctionReturn(0); 277 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 278 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 279 if (yy == zz) { 280 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 281 PLogObjectParent(A,tmp); 282 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 283 } 284 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 285 /* assume if pivots exist then use LU; else Cholesky */ 286 if (mat->pivots) { 287 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 288 } else { 289 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 290 } 291 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 292 if (tmp) { 293 ierr = VecAXPY(&sone,tmp,yy); CHKERRQ(ierr); 294 ierr = VecDestroy(tmp); CHKERRQ(ierr); 295 } else { 296 ierr = VecAXPY(&sone,zz,yy); CHKERRQ(ierr); 297 } 298 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 299 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 300 PLogFlops(mat->n*mat->n - mat->n); 301 PetscFunctionReturn(0); 302 } 303 /* ------------------------------------------------------------------*/ 304 #undef __FUNC__ 305 #define __FUNC__ "MatRelax_SeqDense" 306 int MatRelax_SeqDense(Mat A,Vec bb,double omega,MatSORType flag, 307 double shift,int its,Vec xx) 308 { 309 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 310 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 311 int ierr, m = mat->m, i; 312 #if !defined(USE_PETSC_COMPLEX) 313 int o = 1; 314 #endif 315 316 PetscFunctionBegin; 317 if (flag & SOR_ZERO_INITIAL_GUESS) { 318 /* this is a hack fix, should have another version without the second BLdot */ 319 ierr = VecSet(&zero,xx); CHKERRQ(ierr); 320 } 321 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 322 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 323 while (its--) { 324 if (flag & SOR_FORWARD_SWEEP){ 325 for ( i=0; i<m; i++ ) { 326 #if defined(USE_PETSC_COMPLEX) 327 /* cannot use BLAS dot for complex because compiler/linker is 328 not happy about returning a double complex */ 329 int _i; 330 Scalar sum = b[i]; 331 for ( _i=0; _i<m; _i++ ) { 332 sum -= PetscConj(v[i+_i*m])*x[_i]; 333 } 334 xt = sum; 335 #else 336 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 337 #endif 338 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 339 } 340 } 341 if (flag & SOR_BACKWARD_SWEEP) { 342 for ( i=m-1; i>=0; i-- ) { 343 #if defined(USE_PETSC_COMPLEX) 344 /* cannot use BLAS dot for complex because compiler/linker is 345 not happy about returning a double complex */ 346 int _i; 347 Scalar sum = b[i]; 348 for ( _i=0; _i<m; _i++ ) { 349 sum -= PetscConj(v[i+_i*m])*x[_i]; 350 } 351 xt = sum; 352 #else 353 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 354 #endif 355 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 356 } 357 } 358 } 359 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 360 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 361 PetscFunctionReturn(0); 362 } 363 364 /* -----------------------------------------------------------------*/ 365 #undef __FUNC__ 366 #define __FUNC__ "MatMultTrans_SeqDense" 367 int MatMultTrans_SeqDense(Mat A,Vec xx,Vec yy) 368 { 369 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 370 Scalar *v = mat->v, *x, *y; 371 int ierr,_One=1;Scalar _DOne=1.0, _DZero=0.0; 372 373 PetscFunctionBegin; 374 if (!mat->m || !mat->n) PetscFunctionReturn(0); 375 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 376 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 377 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 378 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 379 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 380 PLogFlops(2*mat->m*mat->n - mat->n); 381 PetscFunctionReturn(0); 382 } 383 384 #undef __FUNC__ 385 #define __FUNC__ "MatMult_SeqDense" 386 int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 387 { 388 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 389 Scalar *v = mat->v, *x, *y; 390 int ierr,_One=1;Scalar _DOne=1.0, _DZero=0.0; 391 392 PetscFunctionBegin; 393 if (!mat->m || !mat->n) PetscFunctionReturn(0); 394 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 395 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 396 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DZero,y,&_One); 397 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 398 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 399 PLogFlops(2*mat->m*mat->n - mat->m); 400 PetscFunctionReturn(0); 401 } 402 403 #undef __FUNC__ 404 #define __FUNC__ "MatMultAdd_SeqDense" 405 int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 406 { 407 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 408 Scalar *v = mat->v, *x, *y; 409 int ierr,_One=1; Scalar _DOne=1.0; 410 411 PetscFunctionBegin; 412 if (!mat->m || !mat->n) PetscFunctionReturn(0); 413 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 414 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 415 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 416 LAgemv_( "N", &(mat->m), &(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 417 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 418 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 419 PLogFlops(2*mat->m*mat->n); 420 PetscFunctionReturn(0); 421 } 422 423 #undef __FUNC__ 424 #define __FUNC__ "MatMultTransAdd_SeqDense" 425 int MatMultTransAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 426 { 427 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 428 Scalar *v = mat->v, *x, *y; 429 int ierr,_One=1; 430 Scalar _DOne=1.0; 431 432 PetscFunctionBegin; 433 if (!mat->m || !mat->n) PetscFunctionReturn(0); 434 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 435 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 436 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 437 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DOne,y,&_One); 438 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 439 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 440 PLogFlops(2*mat->m*mat->n); 441 PetscFunctionReturn(0); 442 } 443 444 /* -----------------------------------------------------------------*/ 445 #undef __FUNC__ 446 #define __FUNC__ "MatGetRow_SeqDense" 447 int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 448 { 449 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 450 Scalar *v; 451 int i; 452 453 PetscFunctionBegin; 454 *ncols = mat->n; 455 if (cols) { 456 *cols = (int *) PetscMalloc((mat->n+1)*sizeof(int)); CHKPTRQ(*cols); 457 for ( i=0; i<mat->n; i++ ) (*cols)[i] = i; 458 } 459 if (vals) { 460 *vals = (Scalar *) PetscMalloc((mat->n+1)*sizeof(Scalar)); CHKPTRQ(*vals); 461 v = mat->v + row; 462 for ( i=0; i<mat->n; i++ ) {(*vals)[i] = *v; v += mat->m;} 463 } 464 PetscFunctionReturn(0); 465 } 466 467 #undef __FUNC__ 468 #define __FUNC__ "MatRestoreRow_SeqDense" 469 int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 470 { 471 if (cols) { PetscFree(*cols); } 472 if (vals) { PetscFree(*vals); } 473 PetscFunctionReturn(0); 474 } 475 /* ----------------------------------------------------------------*/ 476 #undef __FUNC__ 477 #define __FUNC__ "MatSetValues_SeqDense" 478 int MatSetValues_SeqDense(Mat A,int m,int *indexm,int n, 479 int *indexn,Scalar *v,InsertMode addv) 480 { 481 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 482 int i,j; 483 484 PetscFunctionBegin; 485 if (!mat->roworiented) { 486 if (addv == INSERT_VALUES) { 487 for ( j=0; j<n; j++ ) { 488 #if defined(USE_PETSC_BOPT_g) 489 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 490 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 491 #endif 492 for ( i=0; i<m; i++ ) { 493 #if defined(USE_PETSC_BOPT_g) 494 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 495 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 496 #endif 497 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 498 } 499 } 500 } else { 501 for ( j=0; j<n; j++ ) { 502 #if defined(USE_PETSC_BOPT_g) 503 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 504 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 505 #endif 506 for ( i=0; i<m; i++ ) { 507 #if defined(USE_PETSC_BOPT_g) 508 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 509 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 510 #endif 511 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 512 } 513 } 514 } 515 } else { 516 if (addv == INSERT_VALUES) { 517 for ( i=0; i<m; i++ ) { 518 #if defined(USE_PETSC_BOPT_g) 519 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 520 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 521 #endif 522 for ( j=0; j<n; j++ ) { 523 #if defined(USE_PETSC_BOPT_g) 524 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 525 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 526 #endif 527 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 528 } 529 } 530 } else { 531 for ( i=0; i<m; i++ ) { 532 #if defined(USE_PETSC_BOPT_g) 533 if (indexm[i] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 534 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 535 #endif 536 for ( j=0; j<n; j++ ) { 537 #if defined(USE_PETSC_BOPT_g) 538 if (indexn[j] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 539 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 540 #endif 541 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 542 } 543 } 544 } 545 } 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNC__ 550 #define __FUNC__ "MatGetValues_SeqDense" 551 int MatGetValues_SeqDense(Mat A,int m,int *indexm,int n,int *indexn,Scalar *v) 552 { 553 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 554 int i, j; 555 Scalar *vpt = v; 556 557 PetscFunctionBegin; 558 /* row-oriented output */ 559 for ( i=0; i<m; i++ ) { 560 for ( j=0; j<n; j++ ) { 561 *vpt++ = mat->v[indexn[j]*mat->m + indexm[i]]; 562 } 563 } 564 PetscFunctionReturn(0); 565 } 566 567 /* -----------------------------------------------------------------*/ 568 569 #include "sys.h" 570 571 #undef __FUNC__ 572 #define __FUNC__ "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 MPI_Comm_size(comm,&size); 584 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"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,0,"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 PetscFree(w); 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 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 632 633 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 634 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 635 } 636 PetscFunctionReturn(0); 637 } 638 639 #include "pinclude/pviewer.h" 640 #include "sys.h" 641 642 #undef __FUNC__ 643 #define __FUNC__ "MatView_SeqDense_ASCII" 644 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 645 { 646 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 647 int ierr, i, j, format; 648 FILE *fd; 649 char *outputname; 650 Scalar *v; 651 652 PetscFunctionBegin; 653 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 654 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 655 ierr = ViewerGetFormat(viewer,&format); 656 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 657 PetscFunctionReturn(0); /* do nothing for now */ 658 } 659 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 660 for ( i=0; i<a->m; i++ ) { 661 v = a->v + i; 662 fprintf(fd,"row %d:",i); 663 for ( j=0; j<a->n; j++ ) { 664 #if defined(USE_PETSC_COMPLEX) 665 if (PetscReal(*v) != 0.0 && PetscImaginary(*v) != 0.0) fprintf(fd," %d %g + %g i",j,PetscReal(*v),PetscImaginary(*v)); 666 else if (PetscReal(*v)) fprintf(fd," %d %g ",j,PetscReal(*v)); 667 #else 668 if (*v) fprintf(fd," %d %g ",j,*v); 669 #endif 670 v += a->m; 671 } 672 fprintf(fd,"\n"); 673 } 674 } else { 675 #if defined(USE_PETSC_COMPLEX) 676 int allreal = 1; 677 /* determine if matrix has all real values */ 678 v = a->v; 679 for ( i=0; i<a->m*a->n; i++ ) { 680 if (PetscImaginary(v[i])) { allreal = 0; break ;} 681 } 682 #endif 683 for ( i=0; i<a->m; i++ ) { 684 v = a->v + i; 685 for ( j=0; j<a->n; j++ ) { 686 #if defined(USE_PETSC_COMPLEX) 687 if (allreal) { 688 fprintf(fd,"%6.4e ",PetscReal(*v)); v += a->m; 689 } else { 690 fprintf(fd,"%6.4e + %6.4e i ",PetscReal(*v),PetscImaginary(*v)); v += a->m; 691 } 692 #else 693 fprintf(fd,"%6.4e ",*v); v += a->m; 694 #endif 695 } 696 fprintf(fd,"\n"); 697 } 698 } 699 fflush(fd); 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNC__ 704 #define __FUNC__ "MatView_SeqDense_Binary" 705 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 706 { 707 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 708 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 709 int format; 710 Scalar *v, *anonz,*vals; 711 712 PetscFunctionBegin; 713 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 714 715 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 716 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 717 /* store the matrix as a dense matrix */ 718 col_lens = (int *) PetscMalloc( 4*sizeof(int) ); CHKPTRQ(col_lens); 719 col_lens[0] = MAT_COOKIE; 720 col_lens[1] = m; 721 col_lens[2] = n; 722 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 723 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1); CHKERRQ(ierr); 724 PetscFree(col_lens); 725 726 /* write out matrix, by rows */ 727 vals = (Scalar *) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(vals); 728 v = a->v; 729 for ( i=0; i<m; i++ ) { 730 for ( j=0; j<n; j++ ) { 731 vals[i + j*m] = *v++; 732 } 733 } 734 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0); CHKERRQ(ierr); 735 PetscFree(vals); 736 } else { 737 col_lens = (int *) PetscMalloc( (4+nz)*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] = nz; 742 743 /* store lengths of each row and write (including header) to file */ 744 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 745 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1); CHKERRQ(ierr); 746 747 /* Possibly should write in smaller increments, not whole matrix at once? */ 748 /* store column indices (zero start index) */ 749 ict = 0; 750 for ( i=0; i<m; i++ ) { 751 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 752 } 753 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0); CHKERRQ(ierr); 754 PetscFree(col_lens); 755 756 /* store nonzero values */ 757 anonz = (Scalar *) PetscMalloc((nz+1)*sizeof(Scalar)); CHKPTRQ(anonz); 758 ict = 0; 759 for ( i=0; i<m; i++ ) { 760 v = a->v + i; 761 for ( j=0; j<n; j++ ) { 762 anonz[ict++] = *v; v += a->m; 763 } 764 } 765 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0); CHKERRQ(ierr); 766 PetscFree(anonz); 767 } 768 PetscFunctionReturn(0); 769 } 770 771 #undef __FUNC__ 772 #define __FUNC__ "MatView_SeqDense" 773 int MatView_SeqDense(Mat A,Viewer viewer) 774 { 775 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 776 ViewerType vtype; 777 int ierr; 778 779 PetscFunctionBegin; 780 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 781 782 if (vtype == MATLAB_VIEWER) { 783 ierr = ViewerMatlabPutScalar_Private(viewer,a->m,a->n,a->v); CHKERRQ(ierr); 784 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 785 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 786 } else if (vtype == BINARY_FILE_VIEWER) { 787 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 788 } else { 789 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 790 } 791 PetscFunctionReturn(0); 792 } 793 794 #undef __FUNC__ 795 #define __FUNC__ "MatDestroy_SeqDense" 796 int MatDestroy_SeqDense(Mat mat) 797 { 798 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 799 int ierr; 800 801 PetscFunctionBegin; 802 if (--mat->refct > 0) PetscFunctionReturn(0); 803 804 if (mat->mapping) { 805 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 806 } 807 if (mat->bmapping) { 808 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping); CHKERRQ(ierr); 809 } 810 #if defined(USE_PETSC_LOG) 811 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 812 #endif 813 if (l->pivots) PetscFree(l->pivots); 814 if (!l->user_alloc) PetscFree(l->v); 815 PetscFree(l); 816 if (mat->rmap) { 817 ierr = MapDestroy(mat->rmap);CHKERRQ(ierr); 818 } 819 if (mat->cmap) { 820 ierr = MapDestroy(mat->cmap);CHKERRQ(ierr); 821 } 822 PLogObjectDestroy(mat); 823 PetscHeaderDestroy(mat); 824 PetscFunctionReturn(0); 825 } 826 827 #undef __FUNC__ 828 #define __FUNC__ "MatTranspose_SeqDense" 829 int MatTranspose_SeqDense(Mat A,Mat *matout) 830 { 831 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 832 int k, j, m, n; 833 Scalar *v, tmp; 834 835 PetscFunctionBegin; 836 v = mat->v; m = mat->m; n = mat->n; 837 if (matout == PETSC_NULL) { /* in place transpose */ 838 if (m != n) { /* malloc temp to hold transpose */ 839 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 840 for ( j=0; j<m; j++ ) { 841 for ( k=0; k<n; k++ ) { 842 w[k + j*n] = v[j + k*m]; 843 } 844 } 845 PetscMemcpy(v,w,m*n*sizeof(Scalar)); 846 PetscFree(w); 847 } else { 848 for ( j=0; j<m; j++ ) { 849 for ( k=0; k<j; k++ ) { 850 tmp = v[j + k*n]; 851 v[j + k*n] = v[k + j*n]; 852 v[k + j*n] = tmp; 853 } 854 } 855 } 856 } else { /* out-of-place transpose */ 857 int ierr; 858 Mat tmat; 859 Mat_SeqDense *tmatd; 860 Scalar *v2; 861 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 862 tmatd = (Mat_SeqDense *) tmat->data; 863 v = mat->v; v2 = tmatd->v; 864 for ( j=0; j<n; j++ ) { 865 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 866 } 867 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 868 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 869 *matout = tmat; 870 } 871 PetscFunctionReturn(0); 872 } 873 874 #undef __FUNC__ 875 #define __FUNC__ "MatEqual_SeqDense" 876 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 877 { 878 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 879 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 880 int i; 881 Scalar *v1 = mat1->v, *v2 = mat2->v; 882 883 PetscFunctionBegin; 884 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 885 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 886 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 887 for ( i=0; i<mat1->m*mat1->n; i++ ) { 888 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 889 v1++; v2++; 890 } 891 *flg = PETSC_TRUE; 892 PetscFunctionReturn(0); 893 } 894 895 #undef __FUNC__ 896 #define __FUNC__ "MatGetDiagonal_SeqDense" 897 int MatGetDiagonal_SeqDense(Mat A,Vec v) 898 { 899 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 900 int ierr,i, n, len; 901 Scalar *x, zero = 0.0; 902 903 PetscFunctionBegin; 904 ierr = VecSet(&zero,v);CHKERRQ(ierr); 905 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 906 ierr = VecGetArray(v,&x); CHKERRQ(ierr); 907 len = PetscMin(mat->m,mat->n); 908 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 909 for ( i=0; i<len; i++ ) { 910 x[i] = mat->v[i*mat->m + i]; 911 } 912 ierr = VecRestoreArray(v,&x); CHKERRQ(ierr); 913 PetscFunctionReturn(0); 914 } 915 916 #undef __FUNC__ 917 #define __FUNC__ "MatDiagonalScale_SeqDense" 918 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 919 { 920 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 921 Scalar *l,*r,x,*v; 922 int ierr,i,j,m = mat->m, n = mat->n; 923 924 PetscFunctionBegin; 925 if (ll) { 926 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 927 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 928 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 929 for ( i=0; i<m; i++ ) { 930 x = l[i]; 931 v = mat->v + i; 932 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 933 } 934 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 935 PLogFlops(n*m); 936 } 937 if (rr) { 938 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 939 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 940 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 941 for ( i=0; i<n; i++ ) { 942 x = r[i]; 943 v = mat->v + i*m; 944 for ( j=0; j<m; j++ ) { (*v++) *= x;} 945 } 946 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 947 PLogFlops(n*m); 948 } 949 PetscFunctionReturn(0); 950 } 951 952 #undef __FUNC__ 953 #define __FUNC__ "MatNorm_SeqDense" 954 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 955 { 956 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 957 Scalar *v = mat->v; 958 double sum = 0.0; 959 int i, j; 960 961 PetscFunctionBegin; 962 if (type == NORM_FROBENIUS) { 963 for (i=0; i<mat->n*mat->m; i++ ) { 964 #if defined(USE_PETSC_COMPLEX) 965 sum += PetscReal(PetscConj(*v)*(*v)); v++; 966 #else 967 sum += (*v)*(*v); v++; 968 #endif 969 } 970 *norm = sqrt(sum); 971 PLogFlops(2*mat->n*mat->m); 972 } else if (type == NORM_1) { 973 *norm = 0.0; 974 for ( j=0; j<mat->n; j++ ) { 975 sum = 0.0; 976 for ( i=0; i<mat->m; i++ ) { 977 sum += PetscAbsScalar(*v); v++; 978 } 979 if (sum > *norm) *norm = sum; 980 } 981 PLogFlops(mat->n*mat->m); 982 } else if (type == NORM_INFINITY) { 983 *norm = 0.0; 984 for ( j=0; j<mat->m; j++ ) { 985 v = mat->v + j; 986 sum = 0.0; 987 for ( i=0; i<mat->n; i++ ) { 988 sum += PetscAbsScalar(*v); v += mat->m; 989 } 990 if (sum > *norm) *norm = sum; 991 } 992 PLogFlops(mat->n*mat->m); 993 } else { 994 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 995 } 996 PetscFunctionReturn(0); 997 } 998 999 #undef __FUNC__ 1000 #define __FUNC__ "MatSetOption_SeqDense" 1001 int MatSetOption_SeqDense(Mat A,MatOption op) 1002 { 1003 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 1004 1005 PetscFunctionBegin; 1006 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 1007 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1008 else if (op == MAT_ROWS_SORTED || 1009 op == MAT_ROWS_UNSORTED || 1010 op == MAT_COLUMNS_SORTED || 1011 op == MAT_COLUMNS_UNSORTED || 1012 op == MAT_SYMMETRIC || 1013 op == MAT_STRUCTURALLY_SYMMETRIC || 1014 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1015 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1016 op == MAT_NEW_NONZERO_LOCATION_ERROR || 1017 op == MAT_NO_NEW_DIAGONALS || 1018 op == MAT_YES_NEW_DIAGONALS || 1019 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1020 op == MAT_USE_HASH_TABLE) 1021 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1022 else { 1023 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1024 } 1025 PetscFunctionReturn(0); 1026 } 1027 1028 #undef __FUNC__ 1029 #define __FUNC__ "MatZeroEntries_SeqDense" 1030 int MatZeroEntries_SeqDense(Mat A) 1031 { 1032 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1033 1034 PetscFunctionBegin; 1035 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 1036 PetscFunctionReturn(0); 1037 } 1038 1039 #undef __FUNC__ 1040 #define __FUNC__ "MatGetBlockSize_SeqDense" 1041 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1042 { 1043 PetscFunctionBegin; 1044 *bs = 1; 1045 PetscFunctionReturn(0); 1046 } 1047 1048 #undef __FUNC__ 1049 #define __FUNC__ "MatZeroRows_SeqDense" 1050 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1051 { 1052 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1053 int n = l->n, i, j,ierr,N, *rows; 1054 Scalar *slot; 1055 1056 PetscFunctionBegin; 1057 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1058 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1059 for ( i=0; i<N; i++ ) { 1060 slot = l->v + rows[i]; 1061 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1062 } 1063 if (diag) { 1064 for ( i=0; i<N; i++ ) { 1065 slot = l->v + (n+1)*rows[i]; 1066 *slot = *diag; 1067 } 1068 } 1069 ISRestoreIndices(is,&rows); 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #undef __FUNC__ 1074 #define __FUNC__ "MatGetSize_SeqDense" 1075 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1076 { 1077 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1078 1079 PetscFunctionBegin; 1080 *m = mat->m; *n = mat->n; 1081 PetscFunctionReturn(0); 1082 } 1083 1084 #undef __FUNC__ 1085 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1086 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1087 { 1088 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1089 1090 PetscFunctionBegin; 1091 *m = 0; *n = mat->m; 1092 PetscFunctionReturn(0); 1093 } 1094 1095 #undef __FUNC__ 1096 #define __FUNC__ "MatGetArray_SeqDense" 1097 int MatGetArray_SeqDense(Mat A,Scalar **array) 1098 { 1099 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1100 1101 PetscFunctionBegin; 1102 *array = mat->v; 1103 PetscFunctionReturn(0); 1104 } 1105 1106 #undef __FUNC__ 1107 #define __FUNC__ "MatRestoreArray_SeqDense" 1108 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1109 { 1110 PetscFunctionBegin; 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,MatGetSubMatrixCall 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,MatGetSubMatrixCall 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 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 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 .keywords: dense, matrix, LAPACK, BLAS 1288 1289 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1290 @*/ 1291 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1292 { 1293 Mat B; 1294 Mat_SeqDense *b; 1295 int ierr,flg,size; 1296 1297 PetscFunctionBegin; 1298 MPI_Comm_size(comm,&size); 1299 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1300 1301 *A = 0; 1302 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,comm,MatDestroy,MatView); 1303 PLogObjectCreate(B); 1304 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 1305 PetscMemzero(b,sizeof(Mat_SeqDense)); 1306 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1307 B->ops->destroy = MatDestroy_SeqDense; 1308 B->ops->view = MatView_SeqDense; 1309 B->factor = 0; 1310 B->mapping = 0; 1311 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1312 B->data = (void *) b; 1313 1314 b->m = m; B->m = m; B->M = m; 1315 b->n = n; B->n = n; B->N = n; 1316 1317 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1318 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1319 1320 b->pivots = 0; 1321 b->roworiented = 1; 1322 if (data == PETSC_NULL) { 1323 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 1324 PetscMemzero(b->v,m*n*sizeof(Scalar)); 1325 b->user_alloc = 0; 1326 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1327 } else { /* user-allocated storage */ 1328 b->v = data; 1329 b->user_alloc = 1; 1330 } 1331 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 1332 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 1333 1334 *A = B; 1335 PetscFunctionReturn(0); 1336 } 1337 1338 1339 1340 1341 1342 1343