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