1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.171 1999/06/08 22:55:39 balay 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,ierr; 164 165 PetscFunctionBegin; 166 if (mat->pivots) { 167 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 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(PETSC_USE_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(PETSC_USE_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(PETSC_USE_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 int ierr; 471 PetscFunctionBegin; 472 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 473 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 474 PetscFunctionReturn(0); 475 } 476 /* ----------------------------------------------------------------*/ 477 #undef __FUNC__ 478 #define __FUNC__ "MatSetValues_SeqDense" 479 int MatSetValues_SeqDense(Mat A,int m,int *indexm,int n, 480 int *indexn,Scalar *v,InsertMode addv) 481 { 482 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 483 int i,j; 484 485 PetscFunctionBegin; 486 if (!mat->roworiented) { 487 if (addv == INSERT_VALUES) { 488 for ( j=0; j<n; j++ ) { 489 if (indexn[j] < 0) {v += m; continue;} 490 #if defined(PETSC_USE_BOPT_g) 491 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 492 #endif 493 for ( i=0; i<m; i++ ) { 494 if (indexm[i] < 0) {v++; continue;} 495 #if defined(PETSC_USE_BOPT_g) 496 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 497 #endif 498 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 499 } 500 } 501 } else { 502 for ( j=0; j<n; j++ ) { 503 if (indexn[j] < 0) {v += m; continue;} 504 #if defined(PETSC_USE_BOPT_g) 505 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 506 #endif 507 for ( i=0; i<m; i++ ) { 508 if (indexm[i] < 0) {v++; continue;} 509 #if defined(PETSC_USE_BOPT_g) 510 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 511 #endif 512 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 513 } 514 } 515 } 516 } else { 517 if (addv == INSERT_VALUES) { 518 for ( i=0; i<m; i++ ) { 519 if (indexm[i] < 0) { v += n; continue;} 520 #if defined(PETSC_USE_BOPT_g) 521 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 522 #endif 523 for ( j=0; j<n; j++ ) { 524 if (indexn[j] < 0) { v++; continue;} 525 #if defined(PETSC_USE_BOPT_g) 526 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 527 #endif 528 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 529 } 530 } 531 } else { 532 for ( i=0; i<m; i++ ) { 533 if (indexm[i] < 0) { v += n; continue;} 534 #if defined(PETSC_USE_BOPT_g) 535 if (indexm[i] >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 536 #endif 537 for ( j=0; j<n; j++ ) { 538 if (indexn[j] < 0) { v++; continue;} 539 #if defined(PETSC_USE_BOPT_g) 540 if (indexn[j] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 541 #endif 542 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 543 } 544 } 545 } 546 } 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNC__ 551 #define __FUNC__ "MatGetValues_SeqDense" 552 int MatGetValues_SeqDense(Mat A,int m,int *indexm,int n,int *indexn,Scalar *v) 553 { 554 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 555 int i, j; 556 Scalar *vpt = v; 557 558 PetscFunctionBegin; 559 /* row-oriented output */ 560 for ( i=0; i<m; i++ ) { 561 for ( j=0; j<n; j++ ) { 562 *vpt++ = mat->v[indexn[j]*mat->m + indexm[i]]; 563 } 564 } 565 PetscFunctionReturn(0); 566 } 567 568 /* -----------------------------------------------------------------*/ 569 570 #include "sys.h" 571 572 #undef __FUNC__ 573 #define __FUNC__ "MatLoad_SeqDense" 574 int MatLoad_SeqDense(Viewer viewer,MatType type,Mat *A) 575 { 576 Mat_SeqDense *a; 577 Mat B; 578 int *scols, i, j, nz, ierr, fd, header[4], size; 579 int *rowlengths = 0, M, N, *cols; 580 Scalar *vals, *svals, *v,*w; 581 MPI_Comm comm = ((PetscObject)viewer)->comm; 582 583 PetscFunctionBegin; 584 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 585 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 586 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 587 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 588 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Not matrix object"); 589 M = header[1]; N = header[2]; nz = header[3]; 590 591 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 592 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 593 B = *A; 594 a = (Mat_SeqDense *) B->data; 595 v = a->v; 596 /* Allocate some temp space to read in the values and then flip them 597 from row major to column major */ 598 w = (Scalar *) PetscMalloc((M*N+1)*sizeof(Scalar));CHKPTRQ(w); 599 /* read in nonzero values */ 600 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 601 /* now flip the values and store them in the matrix*/ 602 for ( j=0; j<N; j++ ) { 603 for ( i=0; i<M; i++ ) { 604 *v++ =w[i*N+j]; 605 } 606 } 607 ierr = PetscFree(w);CHKERRQ(ierr); 608 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 609 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 610 } else { 611 /* read row lengths */ 612 rowlengths = (int*) PetscMalloc( (M+1)*sizeof(int) );CHKPTRQ(rowlengths); 613 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 614 615 /* create our matrix */ 616 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 617 B = *A; 618 a = (Mat_SeqDense *) B->data; 619 v = a->v; 620 621 /* read column indices and nonzeros */ 622 cols = scols = (int *) PetscMalloc( (nz+1)*sizeof(int) );CHKPTRQ(cols); 623 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 624 vals = svals = (Scalar *) PetscMalloc( (nz+1)*sizeof(Scalar) );CHKPTRQ(vals); 625 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 626 627 /* insert into matrix */ 628 for ( i=0; i<M; i++ ) { 629 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 630 svals += rowlengths[i]; scols += rowlengths[i]; 631 } 632 ierr = PetscFree(vals);CHKERRQ(ierr); 633 ierr = PetscFree(cols);CHKERRQ(ierr); 634 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 635 636 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 637 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 638 } 639 PetscFunctionReturn(0); 640 } 641 642 #include "sys.h" 643 644 #undef __FUNC__ 645 #define __FUNC__ "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 FILE *fd; 651 char *outputname; 652 Scalar *v; 653 654 PetscFunctionBegin; 655 ierr = ViewerASCIIGetPointer(viewer,&fd);CHKERRQ(ierr); 656 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 657 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 658 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 659 PetscFunctionReturn(0); /* do nothing for now */ 660 } 661 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 662 for ( i=0; i<a->m; i++ ) { 663 v = a->v + i; 664 fprintf(fd,"row %d:",i); 665 for ( j=0; j<a->n; j++ ) { 666 #if defined(PETSC_USE_COMPLEX) 667 if (PetscReal(*v) != 0.0 && PetscImaginary(*v) != 0.0) fprintf(fd," %d %g + %g i",j,PetscReal(*v),PetscImaginary(*v)); 668 else if (PetscReal(*v)) fprintf(fd," %d %g ",j,PetscReal(*v)); 669 #else 670 if (*v) fprintf(fd," %d %g ",j,*v); 671 #endif 672 v += a->m; 673 } 674 fprintf(fd,"\n"); 675 } 676 } else { 677 #if defined(PETSC_USE_COMPLEX) 678 int allreal = 1; 679 /* determine if matrix has all real values */ 680 v = a->v; 681 for ( i=0; i<a->m*a->n; i++ ) { 682 if (PetscImaginary(v[i])) { allreal = 0; break ;} 683 } 684 #endif 685 for ( i=0; i<a->m; i++ ) { 686 v = a->v + i; 687 for ( j=0; j<a->n; j++ ) { 688 #if defined(PETSC_USE_COMPLEX) 689 if (allreal) { 690 fprintf(fd,"%6.4e ",PetscReal(*v)); v += a->m; 691 } else { 692 fprintf(fd,"%6.4e + %6.4e i ",PetscReal(*v),PetscImaginary(*v)); v += a->m; 693 } 694 #else 695 fprintf(fd,"%6.4e ",*v); v += a->m; 696 #endif 697 } 698 fprintf(fd,"\n"); 699 } 700 } 701 fflush(fd); 702 PetscFunctionReturn(0); 703 } 704 705 #undef __FUNC__ 706 #define __FUNC__ "MatView_SeqDense_Binary" 707 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 708 { 709 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 710 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 711 int format; 712 Scalar *v, *anonz,*vals; 713 714 PetscFunctionBegin; 715 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 716 717 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 718 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 719 /* store the matrix as a dense matrix */ 720 col_lens = (int *) PetscMalloc( 4*sizeof(int) );CHKPTRQ(col_lens); 721 col_lens[0] = MAT_COOKIE; 722 col_lens[1] = m; 723 col_lens[2] = n; 724 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 725 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 726 ierr = PetscFree(col_lens);CHKERRQ(ierr); 727 728 /* write out matrix, by rows */ 729 vals = (Scalar *) PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(vals); 730 v = a->v; 731 for ( i=0; i<m; i++ ) { 732 for ( j=0; j<n; j++ ) { 733 vals[i + j*m] = *v++; 734 } 735 } 736 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 737 ierr = PetscFree(vals);CHKERRQ(ierr); 738 } else { 739 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) );CHKPTRQ(col_lens); 740 col_lens[0] = MAT_COOKIE; 741 col_lens[1] = m; 742 col_lens[2] = n; 743 col_lens[3] = nz; 744 745 /* store lengths of each row and write (including header) to file */ 746 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 747 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 748 749 /* Possibly should write in smaller increments, not whole matrix at once? */ 750 /* store column indices (zero start index) */ 751 ict = 0; 752 for ( i=0; i<m; i++ ) { 753 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 754 } 755 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 756 ierr = PetscFree(col_lens);CHKERRQ(ierr); 757 758 /* store nonzero values */ 759 anonz = (Scalar *) PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(anonz); 760 ict = 0; 761 for ( i=0; i<m; i++ ) { 762 v = a->v + i; 763 for ( j=0; j<n; j++ ) { 764 anonz[ict++] = *v; v += a->m; 765 } 766 } 767 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 768 ierr = PetscFree(anonz);CHKERRQ(ierr); 769 } 770 PetscFunctionReturn(0); 771 } 772 773 #undef __FUNC__ 774 #define __FUNC__ "MatView_SeqDense" 775 int MatView_SeqDense(Mat A,Viewer viewer) 776 { 777 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 778 ViewerType vtype; 779 int ierr; 780 781 PetscFunctionBegin; 782 ierr = ViewerGetType(viewer,&vtype);CHKERRQ(ierr); 783 784 if (PetscTypeCompare(vtype,SOCKET_VIEWER)) { 785 ierr = ViewerSocketPutScalar_Private(viewer,a->m,a->n,a->v);CHKERRQ(ierr); 786 } else if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 787 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 788 } else if (PetscTypeCompare(vtype,BINARY_VIEWER)) { 789 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 790 } else { 791 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 792 } 793 PetscFunctionReturn(0); 794 } 795 796 #undef __FUNC__ 797 #define __FUNC__ "MatDestroy_SeqDense" 798 int MatDestroy_SeqDense(Mat mat) 799 { 800 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 801 int ierr; 802 803 PetscFunctionBegin; 804 if (--mat->refct > 0) PetscFunctionReturn(0); 805 806 if (mat->mapping) { 807 ierr = ISLocalToGlobalMappingDestroy(mat->mapping);CHKERRQ(ierr); 808 } 809 if (mat->bmapping) { 810 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping);CHKERRQ(ierr); 811 } 812 #if defined(PETSC_USE_LOG) 813 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 814 #endif 815 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 816 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 817 ierr = PetscFree(l);CHKERRQ(ierr); 818 if (mat->rmap) { 819 ierr = MapDestroy(mat->rmap);CHKERRQ(ierr); 820 } 821 if (mat->cmap) { 822 ierr = MapDestroy(mat->cmap);CHKERRQ(ierr); 823 } 824 PLogObjectDestroy(mat); 825 PetscHeaderDestroy(mat); 826 PetscFunctionReturn(0); 827 } 828 829 #undef __FUNC__ 830 #define __FUNC__ "MatTranspose_SeqDense" 831 int MatTranspose_SeqDense(Mat A,Mat *matout) 832 { 833 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 834 int k, j, m, n, ierr; 835 Scalar *v, tmp; 836 837 PetscFunctionBegin; 838 v = mat->v; m = mat->m; n = mat->n; 839 if (matout == PETSC_NULL) { /* in place transpose */ 840 if (m != n) { /* malloc temp to hold transpose */ 841 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 842 for ( j=0; j<m; j++ ) { 843 for ( k=0; k<n; k++ ) { 844 w[k + j*n] = v[j + k*m]; 845 } 846 } 847 ierr = PetscMemcpy(v,w,m*n*sizeof(Scalar));CHKERRQ(ierr); 848 ierr = PetscFree(w);CHKERRQ(ierr); 849 } else { 850 for ( j=0; j<m; j++ ) { 851 for ( k=0; k<j; k++ ) { 852 tmp = v[j + k*n]; 853 v[j + k*n] = v[k + j*n]; 854 v[k + j*n] = tmp; 855 } 856 } 857 } 858 } else { /* out-of-place transpose */ 859 Mat tmat; 860 Mat_SeqDense *tmatd; 861 Scalar *v2; 862 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 863 tmatd = (Mat_SeqDense *) tmat->data; 864 v = mat->v; v2 = tmatd->v; 865 for ( j=0; j<n; j++ ) { 866 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 867 } 868 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 869 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 870 *matout = tmat; 871 } 872 PetscFunctionReturn(0); 873 } 874 875 #undef __FUNC__ 876 #define __FUNC__ "MatEqual_SeqDense" 877 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 878 { 879 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 880 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 881 int i; 882 Scalar *v1 = mat1->v, *v2 = mat2->v; 883 884 PetscFunctionBegin; 885 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 886 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 887 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 888 for ( i=0; i<mat1->m*mat1->n; i++ ) { 889 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 890 v1++; v2++; 891 } 892 *flg = PETSC_TRUE; 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNC__ 897 #define __FUNC__ "MatGetDiagonal_SeqDense" 898 int MatGetDiagonal_SeqDense(Mat A,Vec v) 899 { 900 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 901 int ierr,i, n, len; 902 Scalar *x, zero = 0.0; 903 904 PetscFunctionBegin; 905 ierr = VecSet(&zero,v);CHKERRQ(ierr); 906 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 907 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 908 len = PetscMin(mat->m,mat->n); 909 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 910 for ( i=0; i<len; i++ ) { 911 x[i] = mat->v[i*mat->m + i]; 912 } 913 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNC__ 918 #define __FUNC__ "MatDiagonalScale_SeqDense" 919 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 920 { 921 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 922 Scalar *l,*r,x,*v; 923 int ierr,i,j,m = mat->m, n = mat->n; 924 925 PetscFunctionBegin; 926 if (ll) { 927 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 928 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 929 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 930 for ( i=0; i<m; i++ ) { 931 x = l[i]; 932 v = mat->v + i; 933 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 934 } 935 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 936 PLogFlops(n*m); 937 } 938 if (rr) { 939 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 940 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 941 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 942 for ( i=0; i<n; i++ ) { 943 x = r[i]; 944 v = mat->v + i*m; 945 for ( j=0; j<m; j++ ) { (*v++) *= x;} 946 } 947 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 948 PLogFlops(n*m); 949 } 950 PetscFunctionReturn(0); 951 } 952 953 #undef __FUNC__ 954 #define __FUNC__ "MatNorm_SeqDense" 955 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 956 { 957 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 958 Scalar *v = mat->v; 959 double sum = 0.0; 960 int i, j; 961 962 PetscFunctionBegin; 963 if (type == NORM_FROBENIUS) { 964 for (i=0; i<mat->n*mat->m; i++ ) { 965 #if defined(PETSC_USE_COMPLEX) 966 sum += PetscReal(PetscConj(*v)*(*v)); v++; 967 #else 968 sum += (*v)*(*v); v++; 969 #endif 970 } 971 *norm = sqrt(sum); 972 PLogFlops(2*mat->n*mat->m); 973 } else if (type == NORM_1) { 974 *norm = 0.0; 975 for ( j=0; j<mat->n; j++ ) { 976 sum = 0.0; 977 for ( i=0; i<mat->m; i++ ) { 978 sum += PetscAbsScalar(*v); v++; 979 } 980 if (sum > *norm) *norm = sum; 981 } 982 PLogFlops(mat->n*mat->m); 983 } else if (type == NORM_INFINITY) { 984 *norm = 0.0; 985 for ( j=0; j<mat->m; j++ ) { 986 v = mat->v + j; 987 sum = 0.0; 988 for ( i=0; i<mat->n; i++ ) { 989 sum += PetscAbsScalar(*v); v += mat->m; 990 } 991 if (sum > *norm) *norm = sum; 992 } 993 PLogFlops(mat->n*mat->m); 994 } else { 995 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 996 } 997 PetscFunctionReturn(0); 998 } 999 1000 #undef __FUNC__ 1001 #define __FUNC__ "MatSetOption_SeqDense" 1002 int MatSetOption_SeqDense(Mat A,MatOption op) 1003 { 1004 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 1005 1006 PetscFunctionBegin; 1007 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 1008 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1009 else if (op == MAT_ROWS_SORTED || 1010 op == MAT_ROWS_UNSORTED || 1011 op == MAT_COLUMNS_SORTED || 1012 op == MAT_COLUMNS_UNSORTED || 1013 op == MAT_SYMMETRIC || 1014 op == MAT_STRUCTURALLY_SYMMETRIC || 1015 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1016 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1017 op == MAT_NEW_NONZERO_LOCATION_ERR || 1018 op == MAT_NO_NEW_DIAGONALS || 1019 op == MAT_YES_NEW_DIAGONALS || 1020 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1021 op == MAT_USE_HASH_TABLE) 1022 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1023 else { 1024 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1025 } 1026 PetscFunctionReturn(0); 1027 } 1028 1029 #undef __FUNC__ 1030 #define __FUNC__ "MatZeroEntries_SeqDense" 1031 int MatZeroEntries_SeqDense(Mat A) 1032 { 1033 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1034 int ierr; 1035 1036 PetscFunctionBegin; 1037 ierr = PetscMemzero(l->v,l->m*l->n*sizeof(Scalar));CHKERRQ(ierr); 1038 PetscFunctionReturn(0); 1039 } 1040 1041 #undef __FUNC__ 1042 #define __FUNC__ "MatGetBlockSize_SeqDense" 1043 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1044 { 1045 PetscFunctionBegin; 1046 *bs = 1; 1047 PetscFunctionReturn(0); 1048 } 1049 1050 #undef __FUNC__ 1051 #define __FUNC__ "MatZeroRows_SeqDense" 1052 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1053 { 1054 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1055 int n = l->n, i, j,ierr,N, *rows; 1056 Scalar *slot; 1057 1058 PetscFunctionBegin; 1059 ierr = ISGetSize(is,&N);CHKERRQ(ierr); 1060 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1061 for ( i=0; i<N; i++ ) { 1062 slot = l->v + rows[i]; 1063 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1064 } 1065 if (diag) { 1066 for ( i=0; i<N; i++ ) { 1067 slot = l->v + (n+1)*rows[i]; 1068 *slot = *diag; 1069 } 1070 } 1071 ISRestoreIndices(is,&rows); 1072 PetscFunctionReturn(0); 1073 } 1074 1075 #undef __FUNC__ 1076 #define __FUNC__ "MatGetSize_SeqDense" 1077 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1078 { 1079 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1080 1081 PetscFunctionBegin; 1082 *m = mat->m; *n = mat->n; 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNC__ 1087 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1088 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1089 { 1090 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1091 1092 PetscFunctionBegin; 1093 *m = 0; *n = mat->m; 1094 PetscFunctionReturn(0); 1095 } 1096 1097 #undef __FUNC__ 1098 #define __FUNC__ "MatGetArray_SeqDense" 1099 int MatGetArray_SeqDense(Mat A,Scalar **array) 1100 { 1101 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1102 1103 PetscFunctionBegin; 1104 *array = mat->v; 1105 PetscFunctionReturn(0); 1106 } 1107 1108 #undef __FUNC__ 1109 #define __FUNC__ "MatRestoreArray_SeqDense" 1110 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1111 { 1112 PetscFunctionBegin; 1113 *array = 0; /* user cannot accidently use the array later */ 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNC__ 1118 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1119 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1120 { 1121 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1122 int i, j, ierr, m = mat->m, *irow, *icol, nrows, ncols; 1123 Scalar *av, *bv, *v = mat->v; 1124 Mat newmat; 1125 1126 PetscFunctionBegin; 1127 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1128 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1129 ierr = ISGetSize(isrow,&nrows);CHKERRQ(ierr); 1130 ierr = ISGetSize(iscol,&ncols);CHKERRQ(ierr); 1131 1132 /* Check submatrixcall */ 1133 if (scall == MAT_REUSE_MATRIX) { 1134 int n_cols,n_rows; 1135 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1136 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1137 newmat = *B; 1138 } else { 1139 /* Create and fill new matrix */ 1140 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1141 } 1142 1143 /* Now extract the data pointers and do the copy, column at a time */ 1144 bv = ((Mat_SeqDense*)newmat->data)->v; 1145 1146 for ( i=0; i<ncols; i++ ) { 1147 av = v + m*icol[i]; 1148 for (j=0; j<nrows; j++ ) { 1149 *bv++ = av[irow[j]]; 1150 } 1151 } 1152 1153 /* Assemble the matrices so that the correct flags are set */ 1154 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1155 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1156 1157 /* Free work space */ 1158 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1159 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1160 *B = newmat; 1161 PetscFunctionReturn(0); 1162 } 1163 1164 #undef __FUNC__ 1165 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1166 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B) 1167 { 1168 int ierr,i; 1169 1170 PetscFunctionBegin; 1171 if (scall == MAT_INITIAL_MATRIX) { 1172 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) );CHKPTRQ(*B); 1173 } 1174 1175 for ( i=0; i<n; i++ ) { 1176 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1177 } 1178 PetscFunctionReturn(0); 1179 } 1180 1181 #undef __FUNC__ 1182 #define __FUNC__ "MatCopy_SeqDense" 1183 int MatCopy_SeqDense(Mat A, Mat B,MatStructure str) 1184 { 1185 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1186 int ierr; 1187 1188 PetscFunctionBegin; 1189 if (B->type != MATSEQDENSE) { 1190 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1194 ierr = PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar));CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /* -------------------------------------------------------------------*/ 1199 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1200 MatGetRow_SeqDense, 1201 MatRestoreRow_SeqDense, 1202 MatMult_SeqDense, 1203 MatMultAdd_SeqDense, 1204 MatMultTrans_SeqDense, 1205 MatMultTransAdd_SeqDense, 1206 MatSolve_SeqDense, 1207 MatSolveAdd_SeqDense, 1208 MatSolveTrans_SeqDense, 1209 MatSolveTransAdd_SeqDense, 1210 MatLUFactor_SeqDense, 1211 MatCholeskyFactor_SeqDense, 1212 MatRelax_SeqDense, 1213 MatTranspose_SeqDense, 1214 MatGetInfo_SeqDense, 1215 MatEqual_SeqDense, 1216 MatGetDiagonal_SeqDense, 1217 MatDiagonalScale_SeqDense, 1218 MatNorm_SeqDense, 1219 0, 1220 0, 1221 0, 1222 MatSetOption_SeqDense, 1223 MatZeroEntries_SeqDense, 1224 MatZeroRows_SeqDense, 1225 MatLUFactorSymbolic_SeqDense, 1226 MatLUFactorNumeric_SeqDense, 1227 MatCholeskyFactorSymbolic_SeqDense, 1228 MatCholeskyFactorNumeric_SeqDense, 1229 MatGetSize_SeqDense, 1230 MatGetSize_SeqDense, 1231 MatGetOwnershipRange_SeqDense, 1232 0, 1233 0, 1234 MatGetArray_SeqDense, 1235 MatRestoreArray_SeqDense, 1236 MatDuplicate_SeqDense, 1237 0, 1238 0, 1239 0, 1240 0, 1241 MatAXPY_SeqDense, 1242 MatGetSubMatrices_SeqDense, 1243 0, 1244 MatGetValues_SeqDense, 1245 MatCopy_SeqDense, 1246 0, 1247 MatScale_SeqDense, 1248 0, 1249 0, 1250 0, 1251 MatGetBlockSize_SeqDense, 1252 0, 1253 0, 1254 0, 1255 0, 1256 0, 1257 0, 1258 0, 1259 0, 1260 0, 1261 0, 1262 0, 1263 0, 1264 MatGetMaps_Petsc}; 1265 1266 #undef __FUNC__ 1267 #define __FUNC__ "MatCreateSeqDense" 1268 /*@C 1269 MatCreateSeqDense - Creates a sequential dense matrix that 1270 is stored in column major order (the usual Fortran 77 manner). Many 1271 of the matrix operations use the BLAS and LAPACK routines. 1272 1273 Collective on MPI_Comm 1274 1275 Input Parameters: 1276 + comm - MPI communicator, set to PETSC_COMM_SELF 1277 . m - number of rows 1278 . n - number of columns 1279 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1280 to control all matrix memory allocation. 1281 1282 Output Parameter: 1283 . A - the matrix 1284 1285 Notes: 1286 The data input variable is intended primarily for Fortran programmers 1287 who wish to allocate their own matrix memory space. Most users should 1288 set data=PETSC_NULL. 1289 1290 Level: intermediate 1291 1292 .keywords: dense, matrix, LAPACK, BLAS 1293 1294 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1295 @*/ 1296 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1297 { 1298 Mat B; 1299 Mat_SeqDense *b; 1300 int ierr,flg,size; 1301 1302 PetscFunctionBegin; 1303 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1304 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1305 1306 *A = 0; 1307 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,"Mat",comm,MatDestroy,MatView); 1308 PLogObjectCreate(B); 1309 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense));CHKPTRQ(b); 1310 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1311 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1312 B->ops->destroy = MatDestroy_SeqDense; 1313 B->ops->view = MatView_SeqDense; 1314 B->factor = 0; 1315 B->mapping = 0; 1316 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1317 B->data = (void *) b; 1318 1319 b->m = m; B->m = m; B->M = m; 1320 b->n = n; B->n = n; B->N = n; 1321 1322 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1323 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1324 1325 b->pivots = 0; 1326 b->roworiented = 1; 1327 if (data == PETSC_NULL) { 1328 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(b->v); 1329 ierr = PetscMemzero(b->v,m*n*sizeof(Scalar));CHKERRQ(ierr); 1330 b->user_alloc = 0; 1331 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1332 } else { /* user-allocated storage */ 1333 b->v = data; 1334 b->user_alloc = 1; 1335 } 1336 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 1337 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr);} 1338 1339 *A = B; 1340 PetscFunctionReturn(0); 1341 } 1342 1343 1344 1345 1346 1347 1348