1 #ifndef lint 2 static char vcid[] = "$Id: dense.c,v 1.73 1995/11/03 02:44:58 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Defines the basic matrix operations for sequential dense. 6 */ 7 8 #include "dense.h" 9 #include "pinclude/plapack.h" 10 #include "pinclude/pviewer.h" 11 12 int MatAXPY_SeqDense(Scalar *alpha,Mat X,Mat Y) 13 { 14 Mat_SeqDense *x = (Mat_SeqDense*) X->data,*y = (Mat_SeqDense*) Y->data; 15 int N = x->m*x->n, one = 1; 16 BLaxpy_( &N, alpha, x->v, &one, y->v, &one ); 17 PLogFlops(2*N-1); 18 return 0; 19 } 20 21 static int MatGetInfo_SeqDense(Mat A,MatInfoType flag,int *nz,int *nzalloc,int *mem) 22 { 23 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 24 int i,N = mat->m*mat->n,count = 0; 25 Scalar *v = mat->v; 26 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 27 *nz = count; *nzalloc = N; *mem = (int)A->mem; 28 return 0; 29 } 30 31 /* ---------------------------------------------------------------*/ 32 /* COMMENT: I have chosen to hide column permutation in the pivots, 33 rather than put it in the Mat->col slot.*/ 34 static int MatLUFactor_SeqDense(Mat A,IS row,IS col,double f) 35 { 36 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 37 int info; 38 39 if (!mat->pivots) { 40 mat->pivots = (int *) PetscMalloc(mat->m*sizeof(int));CHKPTRQ(mat->pivots); 41 PLogObjectMemory(A,mat->m*sizeof(int)); 42 } 43 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 44 if (info) SETERRQ(1,"MatLUFactor_SeqDense:Bad LU factorization"); 45 A->factor = FACTOR_LU; 46 PLogFlops((2*mat->n*mat->n*mat->n)/3); 47 return 0; 48 } 49 static int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,double f,Mat *fact) 50 { 51 int ierr; 52 ierr = MatConvert(A,MATSAME,fact); CHKERRQ(ierr); 53 return 0; 54 } 55 static int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 56 { 57 return MatLUFactor(*fact,0,0,1.0); 58 } 59 static int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,double f,Mat *fact) 60 { 61 int ierr; 62 ierr = MatConvert(A,MATSAME,fact); CHKERRQ(ierr); 63 return 0; 64 } 65 static int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 66 { 67 return MatCholeskyFactor(*fact,0,1.0); 68 } 69 static int MatCholeskyFactor_SeqDense(Mat A,IS perm,double f) 70 { 71 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 72 int info; 73 74 if (mat->pivots) { 75 PetscFree(mat->pivots); 76 PLogObjectMemory(A,-mat->m*sizeof(int)); 77 mat->pivots = 0; 78 } 79 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 80 if (info) SETERRQ(1,"MatCholeskyFactor_SeqDense:Bad factorization"); 81 A->factor = FACTOR_CHOLESKY; 82 PLogFlops((mat->n*mat->n*mat->n)/3); 83 return 0; 84 } 85 86 static int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 87 { 88 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 89 int one = 1, info; 90 Scalar *x, *y; 91 92 VecGetArray(xx,&x); VecGetArray(yy,&y); 93 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 94 if (A->factor == FACTOR_LU) { 95 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 96 } 97 else if (A->factor == FACTOR_CHOLESKY){ 98 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 99 } 100 else SETERRQ(1,"MatSolve_SeqDense:Matrix must be factored to solve"); 101 if (info) SETERRQ(1,"MatSolve_SeqDense:Bad solve"); 102 PLogFlops(mat->n*mat->n - mat->n); 103 return 0; 104 } 105 static int MatSolveTrans_SeqDense(Mat A,Vec xx,Vec yy) 106 { 107 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 108 int one = 1, info; 109 Scalar *x, *y; 110 111 VecGetArray(xx,&x); VecGetArray(yy,&y); 112 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 113 /* assume if pivots exist then use LU; else Cholesky */ 114 if (mat->pivots) { 115 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 116 } 117 else { 118 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 119 } 120 if (info) SETERRQ(1,"MatSolveTrans_SeqDense:Bad solve"); 121 PLogFlops(mat->n*mat->n - mat->n); 122 return 0; 123 } 124 static int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 125 { 126 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 127 int one = 1, info,ierr; 128 Scalar *x, *y, sone = 1.0; 129 Vec tmp = 0; 130 131 VecGetArray(xx,&x); VecGetArray(yy,&y); 132 if (yy == zz) { 133 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 134 PLogObjectParent(A,tmp); 135 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 136 } 137 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 138 /* assume if pivots exist then use LU; else Cholesky */ 139 if (mat->pivots) { 140 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 141 } 142 else { 143 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 144 } 145 if (info) SETERRQ(1,"MatSolveAdd_SeqDense:Bad solve"); 146 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 147 else VecAXPY(&sone,zz,yy); 148 PLogFlops(mat->n*mat->n - mat->n); 149 return 0; 150 } 151 152 static int MatSolveTransAdd_SeqDense(Mat A,Vec xx,Vec zz, Vec yy) 153 { 154 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 155 int one = 1, info,ierr; 156 Scalar *x, *y, sone = 1.0; 157 Vec tmp; 158 159 VecGetArray(xx,&x); VecGetArray(yy,&y); 160 if (yy == zz) { 161 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 162 PLogObjectParent(A,tmp); 163 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 164 } 165 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 166 /* assume if pivots exist then use LU; else Cholesky */ 167 if (mat->pivots) { 168 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 169 } 170 else { 171 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 172 } 173 if (info) SETERRQ(1,"MatSolveTransAdd_SeqDense:Bad solve"); 174 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 175 else VecAXPY(&sone,zz,yy); 176 PLogFlops(mat->n*mat->n - mat->n); 177 return 0; 178 } 179 /* ------------------------------------------------------------------*/ 180 static int MatRelax_SeqDense(Mat A,Vec bb,double omega,MatSORType flag, 181 double shift,int its,Vec xx) 182 { 183 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 184 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 185 int o = 1,ierr, m = mat->m, i; 186 187 if (flag & SOR_ZERO_INITIAL_GUESS) { 188 /* this is a hack fix, should have another version without 189 the second BLdot */ 190 ierr = VecSet(&zero,xx); CHKERRQ(ierr); 191 } 192 VecGetArray(xx,&x); VecGetArray(bb,&b); 193 while (its--) { 194 if (flag & SOR_FORWARD_SWEEP){ 195 for ( i=0; i<m; i++ ) { 196 #if defined(PETSC_COMPLEX) 197 /* cannot use BLAS dot for complex because compiler/linker is 198 not happy about returning a double complex */ 199 int _i; 200 Scalar sum = b[i]; 201 for ( _i=0; _i<m; _i++ ) { 202 sum -= conj(v[i+_i*m])*x[_i]; 203 } 204 xt = sum; 205 #else 206 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 207 #endif 208 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 209 } 210 } 211 if (flag & SOR_BACKWARD_SWEEP) { 212 for ( i=m-1; i>=0; i-- ) { 213 #if defined(PETSC_COMPLEX) 214 /* cannot use BLAS dot for complex because compiler/linker is 215 not happy about returning a double complex */ 216 int _i; 217 Scalar sum = b[i]; 218 for ( _i=0; _i<m; _i++ ) { 219 sum -= conj(v[i+_i*m])*x[_i]; 220 } 221 xt = sum; 222 #else 223 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 224 #endif 225 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 226 } 227 } 228 } 229 return 0; 230 } 231 232 /* -----------------------------------------------------------------*/ 233 static int MatMultTrans_SeqDense(Mat A,Vec xx,Vec yy) 234 { 235 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 236 Scalar *v = mat->v, *x, *y; 237 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 238 VecGetArray(xx,&x), VecGetArray(yy,&y); 239 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 240 PLogFlops(2*mat->m*mat->n - mat->n); 241 return 0; 242 } 243 static int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 244 { 245 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 246 Scalar *v = mat->v, *x, *y; 247 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 248 VecGetArray(xx,&x); VecGetArray(yy,&y); 249 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DZero,y,&_One); 250 PLogFlops(2*mat->m*mat->n - mat->m); 251 return 0; 252 } 253 static int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 254 { 255 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 256 Scalar *v = mat->v, *x, *y, *z; 257 int _One=1; Scalar _DOne=1.0; 258 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 259 if (zz != yy) PetscMemcpy(y,z,mat->m*sizeof(Scalar)); 260 LAgemv_( "N", &(mat->m), &(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 261 PLogFlops(2*mat->m*mat->n); 262 return 0; 263 } 264 static int MatMultTransAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 265 { 266 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 267 Scalar *v = mat->v, *x, *y, *z; 268 int _One=1; 269 Scalar _DOne=1.0; 270 VecGetArray(xx,&x); VecGetArray(yy,&y); 271 VecGetArray(zz,&z); 272 if (zz != yy) PetscMemcpy(y,z,mat->m*sizeof(Scalar)); 273 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DOne,y,&_One); 274 PLogFlops(2*mat->m*mat->n); 275 return 0; 276 } 277 278 /* -----------------------------------------------------------------*/ 279 static int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 280 { 281 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 282 Scalar *v; 283 int i; 284 285 *ncols = mat->n; 286 if (cols) { 287 *cols = (int *) PetscMalloc(mat->n*sizeof(int)); CHKPTRQ(*cols); 288 for ( i=0; i<mat->n; i++ ) (*cols)[i] = i; 289 } 290 if (vals) { 291 *vals = (Scalar *) PetscMalloc(mat->n*sizeof(Scalar)); CHKPTRQ(*vals); 292 v = mat->v + row; 293 for ( i=0; i<mat->n; i++ ) {(*vals)[i] = *v; v += mat->m;} 294 } 295 return 0; 296 } 297 static int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 298 { 299 if (cols) { PetscFree(*cols); } 300 if (vals) { PetscFree(*vals); } 301 return 0; 302 } 303 /* ----------------------------------------------------------------*/ 304 static int MatInsert_SeqDense(Mat A,int m,int *indexm,int n, 305 int *indexn,Scalar *v,InsertMode addv) 306 { 307 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 308 int i,j; 309 310 if (!mat->roworiented) { 311 if (addv == INSERT_VALUES) { 312 for ( j=0; j<n; j++ ) { 313 for ( i=0; i<m; i++ ) { 314 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 315 } 316 } 317 } 318 else { 319 for ( j=0; j<n; j++ ) { 320 for ( i=0; i<m; i++ ) { 321 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 322 } 323 } 324 } 325 } 326 else { 327 if (addv == INSERT_VALUES) { 328 for ( i=0; i<m; i++ ) { 329 for ( j=0; j<n; j++ ) { 330 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 331 } 332 } 333 } 334 else { 335 for ( i=0; i<m; i++ ) { 336 for ( j=0; j<n; j++ ) { 337 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 338 } 339 } 340 } 341 } 342 return 0; 343 } 344 345 /* -----------------------------------------------------------------*/ 346 static int MatCopyPrivate_SeqDense(Mat A,Mat *newmat,int cpvalues) 347 { 348 Mat_SeqDense *mat = (Mat_SeqDense *) A->data, *l; 349 int ierr; 350 Mat newi; 351 352 ierr = MatCreateSeqDense(A->comm,mat->m,mat->n,&newi); CHKERRQ(ierr); 353 l = (Mat_SeqDense *) newi->data; 354 if (cpvalues == COPY_VALUES) { 355 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 356 } 357 *newmat = newi; 358 return 0; 359 } 360 361 #include "sysio.h" 362 363 int MatLoad_SeqDense(Viewer bview,MatType type,Mat *A) 364 { 365 Mat_SeqDense *a; 366 Mat B; 367 int *scols, i, j, nz, ierr, fd, header[4], size; 368 int *rowlengths = 0, M, N, *cols; 369 Scalar *vals, *svals, *v; 370 PetscObject vobj = (PetscObject) bview; 371 MPI_Comm comm = vobj->comm; 372 373 MPI_Comm_size(comm,&size); 374 if (size > 1) SETERRQ(1,"MatLoad_SeqDense: view must have one processor"); 375 ierr = ViewerFileGetDescriptor_Private(bview,&fd); CHKERRQ(ierr); 376 ierr = SYRead(fd,header,4,SYINT); CHKERRQ(ierr); 377 if (header[0] != MAT_COOKIE) SETERRQ(1,"MatLoad_SeqDense:Not matrix object"); 378 M = header[1]; N = header[2]; nz = header[3]; 379 380 /* read row lengths */ 381 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 382 ierr = SYRead(fd,rowlengths,M,SYINT); CHKERRQ(ierr); 383 384 /* create our matrix */ 385 ierr = MatCreateSeqDense(comm,M,N,A); CHKERRQ(ierr); 386 B = *A; 387 a = (Mat_SeqDense *) B->data; 388 v = a->v; 389 390 /* read column indices and nonzeros */ 391 cols = scols = (int *) PetscMalloc( nz*sizeof(int) ); CHKPTRQ(cols); 392 ierr = SYRead(fd,cols,nz,SYINT); CHKERRQ(ierr); 393 vals = svals = (Scalar *) PetscMalloc( nz*sizeof(Scalar) ); CHKPTRQ(vals); 394 ierr = SYRead(fd,vals,nz,SYSCALAR); CHKERRQ(ierr); 395 396 /* insert into matrix */ 397 for ( i=0; i<M; i++ ) { 398 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 399 svals += rowlengths[i]; scols += rowlengths[i]; 400 } 401 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 402 403 ierr = MatAssemblyBegin(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 404 ierr = MatAssemblyEnd(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 405 return 0; 406 } 407 408 #include "pinclude/pviewer.h" 409 #include "sysio.h" 410 411 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 412 { 413 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 414 int ierr, i, j, format; 415 FILE *fd; 416 char *outputname; 417 Scalar *v; 418 419 ierr = ViewerFileGetPointer_Private(viewer,&fd); CHKERRQ(ierr); 420 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 421 ierr = ViewerFileGetFormat_Private(viewer,&format); 422 if (format == FILE_FORMAT_INFO) { 423 ; /* do nothing for now */ 424 } 425 else { 426 for ( i=0; i<a->m; i++ ) { 427 v = a->v + i; 428 for ( j=0; j<a->n; j++ ) { 429 #if defined(PETSC_COMPLEX) 430 fprintf(fd,"%6.4e + %6.4e i ",real(*v),imag(*v)); v += a->m; 431 #else 432 fprintf(fd,"%6.4e ",*v); v += a->m; 433 #endif 434 } 435 fprintf(fd,"\n"); 436 } 437 } 438 fflush(fd); 439 return 0; 440 } 441 442 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 443 { 444 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 445 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 446 Scalar *v, *anonz; 447 448 ierr = ViewerFileGetDescriptor_Private(viewer,&fd); CHKERRQ(ierr); 449 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) ); CHKPTRQ(col_lens); 450 col_lens[0] = MAT_COOKIE; 451 col_lens[1] = m; 452 col_lens[2] = n; 453 col_lens[3] = nz; 454 455 /* store lengths of each row and write (including header) to file */ 456 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 457 ierr = SYWrite(fd,col_lens,4+m,SYINT,1); CHKERRQ(ierr); 458 459 /* Possibly should write in smaller increments, not whole matrix at once? */ 460 /* store column indices (zero start index) */ 461 ict = 0; 462 for ( i=0; i<m; i++ ) { 463 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 464 } 465 ierr = SYWrite(fd,col_lens,nz,SYINT,0); CHKERRQ(ierr); 466 PetscFree(col_lens); 467 468 /* store nonzero values */ 469 anonz = (Scalar *) PetscMalloc((nz)*sizeof(Scalar)); CHKPTRQ(anonz); 470 ict = 0; 471 for ( i=0; i<m; i++ ) { 472 v = a->v + i; 473 for ( j=0; j<n; j++ ) { 474 anonz[ict++] = *v; v += a->m; 475 } 476 } 477 ierr = SYWrite(fd,anonz,nz,SYSCALAR,0); CHKERRQ(ierr); 478 PetscFree(anonz); 479 return 0; 480 } 481 482 static int MatView_SeqDense(PetscObject obj,Viewer viewer) 483 { 484 Mat A = (Mat) obj; 485 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 486 PetscObject vobj = (PetscObject) viewer; 487 488 if (!viewer) { 489 viewer = STDOUT_VIEWER_SELF; vobj = (PetscObject) viewer; 490 } 491 if (vobj->cookie == VIEWER_COOKIE) { 492 if (vobj->type == MATLAB_VIEWER) { 493 return ViewerMatlabPutArray_Private(viewer,a->m,a->n,a->v); 494 } 495 else if (vobj->type == ASCII_FILE_VIEWER || vobj->type == ASCII_FILES_VIEWER) { 496 return MatView_SeqDense_ASCII(A,viewer); 497 } 498 else if (vobj->type == BINARY_FILE_VIEWER) { 499 return MatView_SeqDense_Binary(A,viewer); 500 } 501 } 502 return 0; 503 } 504 505 static int MatDestroy_SeqDense(PetscObject obj) 506 { 507 Mat mat = (Mat) obj; 508 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 509 #if defined(PETSC_LOG) 510 PLogObjectState(obj,"Rows %d Cols %d",l->m,l->n); 511 #endif 512 if (l->pivots) PetscFree(l->pivots); 513 PetscFree(l); 514 PLogObjectDestroy(mat); 515 PetscHeaderDestroy(mat); 516 return 0; 517 } 518 519 static int MatTranspose_SeqDense(Mat A,Mat *matout) 520 { 521 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 522 int k, j, m, n; 523 Scalar *v, tmp; 524 525 v = mat->v; m = mat->m; n = mat->n; 526 if (!matout) { /* in place transpose */ 527 if (m != n) SETERRQ(1,"MatTranspose_SeqDense:Supports square matrix only in-place"); 528 for ( j=0; j<m; j++ ) { 529 for ( k=0; k<j; k++ ) { 530 tmp = v[j + k*n]; 531 v[j + k*n] = v[k + j*n]; 532 v[k + j*n] = tmp; 533 } 534 } 535 } 536 else { /* out-of-place transpose */ 537 int ierr; 538 Mat tmat; 539 Mat_SeqDense *tmatd; 540 Scalar *v2; 541 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,&tmat); CHKERRQ(ierr); 542 tmatd = (Mat_SeqDense *) tmat->data; 543 v = mat->v; v2 = tmatd->v; 544 for ( j=0; j<n; j++ ) { 545 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 546 } 547 ierr = MatAssemblyBegin(tmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 548 ierr = MatAssemblyEnd(tmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 549 *matout = tmat; 550 } 551 return 0; 552 } 553 554 static int MatEqual_SeqDense(Mat A1,Mat A2) 555 { 556 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 557 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 558 int i; 559 Scalar *v1 = mat1->v, *v2 = mat2->v; 560 if (mat1->m != mat2->m) return 0; 561 if (mat1->n != mat2->n) return 0; 562 for ( i=0; i<mat1->m*mat1->n; i++ ) { 563 if (*v1 != *v2) return 0; 564 v1++; v2++; 565 } 566 return 1; 567 } 568 569 static int MatGetDiagonal_SeqDense(Mat A,Vec v) 570 { 571 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 572 int i, n; 573 Scalar *x; 574 VecGetArray(v,&x); VecGetSize(v,&n); 575 if (n != mat->m) SETERRQ(1,"MatGetDiagonal_SeqDense:Nonconforming mat and vec"); 576 for ( i=0; i<mat->m; i++ ) { 577 x[i] = mat->v[i*mat->m + i]; 578 } 579 return 0; 580 } 581 582 static int MatScale_SeqDense(Mat A,Vec ll,Vec rr) 583 { 584 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 585 Scalar *l,*r,x,*v; 586 int i,j,m = mat->m, n = mat->n; 587 588 if (ll) { 589 VecGetArray(ll,&l); VecGetSize(ll,&m); 590 if (m != mat->m) SETERRQ(1,"MatScale_SeqDense:Left scaling vec wrong size"); 591 PLogFlops(n*m); 592 for ( i=0; i<m; i++ ) { 593 x = l[i]; 594 v = mat->v + i; 595 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 596 } 597 } 598 if (rr) { 599 VecGetArray(rr,&r); VecGetSize(rr,&n); 600 if (n != mat->n) SETERRQ(1,"MatScale_SeqDense:Right scaling vec wrong size"); 601 PLogFlops(n*m); 602 for ( i=0; i<n; i++ ) { 603 x = r[i]; 604 v = mat->v + i*m; 605 for ( j=0; j<m; j++ ) { (*v++) *= x;} 606 } 607 } 608 return 0; 609 } 610 611 static int MatNorm_SeqDense(Mat A,NormType type,double *norm) 612 { 613 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 614 Scalar *v = mat->v; 615 double sum = 0.0; 616 int i, j; 617 618 if (type == NORM_FROBENIUS) { 619 for (i=0; i<mat->n*mat->m; i++ ) { 620 #if defined(PETSC_COMPLEX) 621 sum += real(conj(*v)*(*v)); v++; 622 #else 623 sum += (*v)*(*v); v++; 624 #endif 625 } 626 *norm = sqrt(sum); 627 PLogFlops(2*mat->n*mat->m); 628 } 629 else if (type == NORM_1) { 630 *norm = 0.0; 631 for ( j=0; j<mat->n; j++ ) { 632 sum = 0.0; 633 for ( i=0; i<mat->m; i++ ) { 634 sum += PetscAbsScalar(*v); v++; 635 } 636 if (sum > *norm) *norm = sum; 637 } 638 PLogFlops(mat->n*mat->m); 639 } 640 else if (type == NORM_INFINITY) { 641 *norm = 0.0; 642 for ( j=0; j<mat->m; j++ ) { 643 v = mat->v + j; 644 sum = 0.0; 645 for ( i=0; i<mat->n; i++ ) { 646 sum += PetscAbsScalar(*v); v += mat->m; 647 } 648 if (sum > *norm) *norm = sum; 649 } 650 PLogFlops(mat->n*mat->m); 651 } 652 else { 653 SETERRQ(1,"MatNorm_SeqDense:No two norm"); 654 } 655 return 0; 656 } 657 658 static int MatSetOption_SeqDense(Mat A,MatOption op) 659 { 660 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 661 662 if (op == ROW_ORIENTED) aij->roworiented = 1; 663 else if (op == COLUMN_ORIENTED) aij->roworiented = 0; 664 else if (op == ROWS_SORTED || 665 op == COLUMNS_SORTED || 666 op == SYMMETRIC_MATRIX || 667 op == STRUCTURALLY_SYMMETRIC_MATRIX || 668 op == NO_NEW_NONZERO_LOCATIONS || 669 op == YES_NEW_NONZERO_LOCATIONS || 670 op == NO_NEW_DIAGONALS || 671 op == YES_NEW_DIAGONALS) 672 PLogInfo((PetscObject)A,"Info:MatSetOption_SeqDense:Option ignored\n"); 673 else 674 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqDense:unknown option");} 675 return 0; 676 } 677 678 static int MatZeroEntries_SeqDense(Mat A) 679 { 680 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 681 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 682 return 0; 683 } 684 685 static int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 686 { 687 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 688 int n = l->n, i, j,ierr,N, *rows; 689 Scalar *slot; 690 691 ierr = ISGetLocalSize(is,&N); CHKERRQ(ierr); 692 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 693 for ( i=0; i<N; i++ ) { 694 slot = l->v + rows[i]; 695 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 696 } 697 if (diag) { 698 for ( i=0; i<N; i++ ) { 699 slot = l->v + (n+1)*rows[i]; 700 *slot = *diag; 701 } 702 } 703 ISRestoreIndices(is,&rows); 704 return 0; 705 } 706 707 static int MatGetSize_SeqDense(Mat A,int *m,int *n) 708 { 709 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 710 *m = mat->m; *n = mat->n; 711 return 0; 712 } 713 714 static int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 715 { 716 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 717 *m = 0; *n = mat->m; 718 return 0; 719 } 720 721 static int MatGetArray_SeqDense(Mat A,Scalar **array) 722 { 723 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 724 *array = mat->v; 725 return 0; 726 } 727 728 729 static int MatGetSubMatrixInPlace_SeqDense(Mat A,IS isrow,IS iscol) 730 { 731 SETERRQ(1,"MatGetSubMatrixInPlace_SeqDense:not done"); 732 } 733 734 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall, 735 Mat *submat) 736 { 737 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 738 int nznew, *smap, i, j, ierr, oldcols = mat->n; 739 int *irow, *icol, nrows, ncols, *cwork; 740 Scalar *vwork, *val; 741 Mat newmat; 742 743 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 744 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 745 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 746 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 747 748 smap = (int *) PetscMalloc(oldcols*sizeof(int)); CHKPTRQ(smap); 749 cwork = (int *) PetscMalloc(ncols*sizeof(int)); CHKPTRQ(cwork); 750 vwork = (Scalar *) PetscMalloc(ncols*sizeof(Scalar)); CHKPTRQ(vwork); 751 PetscMemzero((char*)smap,oldcols*sizeof(int)); 752 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 753 754 /* Create and fill new matrix */ 755 ierr = MatCreateSeqDense(A->comm,nrows,ncols,&newmat);CHKERRQ(ierr); 756 for (i=0; i<nrows; i++) { 757 nznew = 0; 758 val = mat->v + irow[i]; 759 for (j=0; j<oldcols; j++) { 760 if (smap[j]) { 761 cwork[nznew] = smap[j] - 1; 762 vwork[nznew++] = val[j * mat->m]; 763 } 764 } 765 ierr = MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES); 766 CHKERRQ(ierr); 767 } 768 ierr = MatAssemblyBegin(newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 769 ierr = MatAssemblyEnd(newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 770 771 /* Free work space */ 772 PetscFree(smap); PetscFree(cwork); PetscFree(vwork); 773 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 774 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 775 *submat = newmat; 776 return 0; 777 } 778 779 /* -------------------------------------------------------------------*/ 780 static struct _MatOps MatOps = {MatInsert_SeqDense, 781 MatGetRow_SeqDense, MatRestoreRow_SeqDense, 782 MatMult_SeqDense, MatMultAdd_SeqDense, 783 MatMultTrans_SeqDense, MatMultTransAdd_SeqDense, 784 MatSolve_SeqDense,MatSolveAdd_SeqDense, 785 MatSolveTrans_SeqDense,MatSolveTransAdd_SeqDense, 786 MatLUFactor_SeqDense,MatCholeskyFactor_SeqDense, 787 MatRelax_SeqDense, 788 MatTranspose_SeqDense, 789 MatGetInfo_SeqDense,MatEqual_SeqDense, 790 MatGetDiagonal_SeqDense,MatScale_SeqDense,MatNorm_SeqDense, 791 0,0, 792 0, MatSetOption_SeqDense,MatZeroEntries_SeqDense,MatZeroRows_SeqDense,0, 793 MatLUFactorSymbolic_SeqDense,MatLUFactorNumeric_SeqDense, 794 MatCholeskyFactorSymbolic_SeqDense,MatCholeskyFactorNumeric_SeqDense, 795 MatGetSize_SeqDense,MatGetSize_SeqDense,MatGetOwnershipRange_SeqDense, 796 0,0,MatGetArray_SeqDense,0,0, 797 MatGetSubMatrix_SeqDense,MatGetSubMatrixInPlace_SeqDense, 798 MatCopyPrivate_SeqDense,0,0,0,0, 799 MatAXPY_SeqDense}; 800 801 /*@C 802 MatCreateSeqDense - Creates a sequential dense matrix that 803 is stored in column major order (the usual Fortran 77 manner). Many 804 of the matrix operations use the BLAS and LAPACK routines. 805 806 Input Parameters: 807 . comm - MPI communicator, set to MPI_COMM_SELF 808 . m - number of rows 809 . n - number of column 810 811 Output Parameter: 812 . newmat - the matrix 813 814 .keywords: dense, matrix, LAPACK, BLAS 815 816 .seealso: MatCreate(), MatSetValues() 817 @*/ 818 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Mat *newmat) 819 { 820 int size = sizeof(Mat_SeqDense) + m*n*sizeof(Scalar); 821 Mat mat; 822 Mat_SeqDense *l; 823 824 *newmat = 0; 825 PetscHeaderCreate(mat,_Mat,MAT_COOKIE,MATSEQDENSE,comm); 826 PLogObjectCreate(mat); 827 l = (Mat_SeqDense *) PetscMalloc(size); CHKPTRQ(l); 828 PetscMemcpy(&mat->ops,&MatOps,sizeof(struct _MatOps)); 829 mat->destroy = MatDestroy_SeqDense; 830 mat->view = MatView_SeqDense; 831 mat->data = (void *) l; 832 mat->factor = 0; 833 PLogObjectMemory(mat,sizeof(struct _Mat) + size); 834 835 l->m = m; 836 l->n = n; 837 l->v = (Scalar *) (l + 1); 838 l->pivots = 0; 839 l->roworiented = 1; 840 841 PetscMemzero(l->v,m*n*sizeof(Scalar)); 842 *newmat = mat; 843 return 0; 844 } 845 846 int MatCreate_SeqDense(Mat A,Mat *newmat) 847 { 848 Mat_SeqDense *m = (Mat_SeqDense *) A->data; 849 return MatCreateSeqDense(A->comm,m->m,m->n,newmat); 850 } 851