1 2 3 #ifndef lint 4 static char vcid[] = "$Id: aij.c,v 1.199 1996/12/02 20:20:36 balay Exp balay $"; 5 #endif 6 7 /* 8 B Defines the basic matrix operations for the AIJ (compressed row) 9 matrix storage format. 10 */ 11 #include "src/mat/impls/aij/seq/aij.h" 12 #include "src/vec/vecimpl.h" 13 #include "src/inline/spops.h" 14 #include "petsc.h" 15 #include "src/inline/bitarray.h" 16 17 /* 18 Basic AIJ format ILU based on drop tolerance 19 */ 20 #undef __FUNCTION__ 21 #define __FUNCTION__ "MatILUDTFactor_SeqAIJ" 22 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 23 { 24 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 25 int ierr = 1; 26 27 SETERRQ(ierr,"MatILUDTFactor_SeqAIJ:Not implemented"); 28 } 29 30 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 31 32 #undef __FUNCTION__ 33 #define __FUNCTION__ "MatGetRowIJ_SeqAIJ" 34 static int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 35 PetscTruth *done) 36 { 37 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 38 int ierr,i,ishift; 39 40 *m = A->m; 41 if (!ia) return 0; 42 ishift = a->indexshift; 43 if (symmetric) { 44 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 45 } else if (oshift == 0 && ishift == -1) { 46 int nz = a->i[a->m]; 47 /* malloc space and subtract 1 from i and j indices */ 48 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 49 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 50 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 51 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 52 } else if (oshift == 1 && ishift == 0) { 53 int nz = a->i[a->m] + 1; 54 /* malloc space and add 1 to i and j indices */ 55 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 56 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 57 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 58 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 59 } else { 60 *ia = a->i; *ja = a->j; 61 } 62 63 return 0; 64 } 65 66 #undef __FUNCTION__ 67 #define __FUNCTION__ "MatRestoreRowIJ_SeqAIJ" 68 static int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 69 PetscTruth *done) 70 { 71 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 72 int ishift = a->indexshift; 73 74 if (!ia) return 0; 75 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 76 PetscFree(*ia); 77 PetscFree(*ja); 78 } 79 return 0; 80 } 81 82 #undef __FUNCTION__ 83 #define __FUNCTION__ "MatGetColumnIJ_SeqAIJ" 84 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 85 PetscTruth *done) 86 { 87 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 88 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 89 int nz = a->i[m]+ishift,row,*jj,mr,col; 90 91 *nn = A->n; 92 if (!ia) return 0; 93 if (symmetric) { 94 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 95 } else { 96 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 97 PetscMemzero(collengths,n*sizeof(int)); 98 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 99 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 100 jj = a->j; 101 for ( i=0; i<nz; i++ ) { 102 collengths[jj[i] + ishift]++; 103 } 104 cia[0] = oshift; 105 for ( i=0; i<n; i++) { 106 cia[i+1] = cia[i] + collengths[i]; 107 } 108 PetscMemzero(collengths,n*sizeof(int)); 109 jj = a->j; 110 for ( row=0; row<m; row++ ) { 111 mr = a->i[row+1] - a->i[row]; 112 for ( i=0; i<mr; i++ ) { 113 col = *jj++ + ishift; 114 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 115 } 116 } 117 PetscFree(collengths); 118 *ia = cia; *ja = cja; 119 } 120 121 return 0; 122 } 123 124 #undef __FUNCTION__ 125 #define __FUNCTION__ "MatRestoreColumnIJ_SeqAIJ" 126 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 127 int **ja,PetscTruth *done) 128 { 129 if (!ia) return 0; 130 131 PetscFree(*ia); 132 PetscFree(*ja); 133 134 return 0; 135 } 136 137 #define CHUNKSIZE 15 138 139 #undef __FUNCTION__ 140 #define __FUNCTION__ "MatSetValues_SeqAIJ" 141 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 142 { 143 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 144 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 145 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 146 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 147 Scalar *ap,value, *aa = a->a; 148 149 for ( k=0; k<m; k++ ) { /* loop over added rows */ 150 row = im[k]; 151 #if defined(PETSC_BOPT_g) 152 if (row < 0) SETERRQ(1,"MatSetValues_SeqAIJ:Negative row"); 153 if (row >= a->m) SETERRQ(1,"MatSetValues_SeqAIJ:Row too large"); 154 #endif 155 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 156 rmax = imax[row]; nrow = ailen[row]; 157 low = 0; 158 for ( l=0; l<n; l++ ) { /* loop over added columns */ 159 #if defined(PETSC_BOPT_g) 160 if (in[l] < 0) SETERRQ(1,"MatSetValues_SeqAIJ:Negative column"); 161 if (in[l] >= a->n) SETERRQ(1,"MatSetValues_SeqAIJ:Column too large"); 162 #endif 163 col = in[l] - shift; 164 if (roworiented) { 165 value = *v++; 166 } 167 else { 168 value = v[k + l*m]; 169 } 170 if (!sorted) low = 0; high = nrow; 171 while (high-low > 5) { 172 t = (low+high)/2; 173 if (rp[t] > col) high = t; 174 else low = t; 175 } 176 for ( i=low; i<high; i++ ) { 177 if (rp[i] > col) break; 178 if (rp[i] == col) { 179 if (is == ADD_VALUES) ap[i] += value; 180 else ap[i] = value; 181 goto noinsert; 182 } 183 } 184 if (nonew) goto noinsert; 185 if (nrow >= rmax) { 186 /* there is no extra room in row, therefore enlarge */ 187 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 188 Scalar *new_a; 189 190 /* malloc new storage space */ 191 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 192 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 193 new_j = (int *) (new_a + new_nz); 194 new_i = new_j + new_nz; 195 196 /* copy over old data into new slots */ 197 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 198 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 199 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 200 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 201 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 202 len*sizeof(int)); 203 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 204 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 205 len*sizeof(Scalar)); 206 /* free up old matrix storage */ 207 PetscFree(a->a); 208 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 209 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 210 a->singlemalloc = 1; 211 212 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 213 rmax = imax[row] = imax[row] + CHUNKSIZE; 214 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 215 a->maxnz += CHUNKSIZE; 216 a->reallocs++; 217 } 218 N = nrow++ - 1; a->nz++; 219 /* shift up all the later entries in this row */ 220 for ( ii=N; ii>=i; ii-- ) { 221 rp[ii+1] = rp[ii]; 222 ap[ii+1] = ap[ii]; 223 } 224 rp[i] = col; 225 ap[i] = value; 226 noinsert:; 227 low = i + 1; 228 } 229 ailen[row] = nrow; 230 } 231 return 0; 232 } 233 234 #undef __FUNCTION__ 235 #define __FUNCTION__ "MatGetValues_SeqAIJ" 236 static int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 237 { 238 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 239 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 240 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 241 Scalar *ap, *aa = a->a, zero = 0.0; 242 243 for ( k=0; k<m; k++ ) { /* loop over rows */ 244 row = im[k]; 245 if (row < 0) SETERRQ(1,"MatGetValues_SeqAIJ:Negative row"); 246 if (row >= a->m) SETERRQ(1,"MatGetValues_SeqAIJ:Row too large"); 247 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 248 nrow = ailen[row]; 249 for ( l=0; l<n; l++ ) { /* loop over columns */ 250 if (in[l] < 0) SETERRQ(1,"MatGetValues_SeqAIJ:Negative column"); 251 if (in[l] >= a->n) SETERRQ(1,"MatGetValues_SeqAIJ:Column too large"); 252 col = in[l] - shift; 253 high = nrow; low = 0; /* assume unsorted */ 254 while (high-low > 5) { 255 t = (low+high)/2; 256 if (rp[t] > col) high = t; 257 else low = t; 258 } 259 for ( i=low; i<high; i++ ) { 260 if (rp[i] > col) break; 261 if (rp[i] == col) { 262 *v++ = ap[i]; 263 goto finished; 264 } 265 } 266 *v++ = zero; 267 finished:; 268 } 269 } 270 return 0; 271 } 272 273 #include "draw.h" 274 #include "pinclude/pviewer.h" 275 #include "sys.h" 276 277 #undef __FUNCTION__ 278 #define __FUNCTION__ "MatView_SeqAIJ_Binary" 279 static int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 280 { 281 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 282 int i, fd, *col_lens, ierr; 283 284 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 285 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 286 col_lens[0] = MAT_COOKIE; 287 col_lens[1] = a->m; 288 col_lens[2] = a->n; 289 col_lens[3] = a->nz; 290 291 /* store lengths of each row and write (including header) to file */ 292 for ( i=0; i<a->m; i++ ) { 293 col_lens[4+i] = a->i[i+1] - a->i[i]; 294 } 295 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,BINARY_INT,1); CHKERRQ(ierr); 296 PetscFree(col_lens); 297 298 /* store column indices (zero start index) */ 299 if (a->indexshift) { 300 for ( i=0; i<a->nz; i++ ) a->j[i]--; 301 } 302 ierr = PetscBinaryWrite(fd,a->j,a->nz,BINARY_INT,0); CHKERRQ(ierr); 303 if (a->indexshift) { 304 for ( i=0; i<a->nz; i++ ) a->j[i]++; 305 } 306 307 /* store nonzero values */ 308 ierr = PetscBinaryWrite(fd,a->a,a->nz,BINARY_SCALAR,0); CHKERRQ(ierr); 309 return 0; 310 } 311 312 #undef __FUNCTION__ 313 #define __FUNCTION__ "MatView_SeqAIJ_ASCII" 314 static int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 315 { 316 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 317 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 318 FILE *fd; 319 char *outputname; 320 321 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 322 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 323 ierr = ViewerGetFormat(viewer,&format); 324 if (format == VIEWER_FORMAT_ASCII_INFO) { 325 return 0; 326 } 327 else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 328 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 329 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 330 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 331 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 332 a->inode.node_count,a->inode.limit); 333 } 334 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 335 fprintf(fd,"%% Size = %d %d \n",m,a->n); 336 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 337 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz); 338 fprintf(fd,"zzz = [\n"); 339 340 for (i=0; i<m; i++) { 341 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 342 #if defined(PETSC_COMPLEX) 343 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,real(a->a[j]), 344 imag(a->a[j])); 345 #else 346 fprintf(fd,"%d %d %18.16e\n", i+1, a->j[j]+!shift, a->a[j]); 347 #endif 348 } 349 } 350 fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 351 } 352 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 353 for ( i=0; i<m; i++ ) { 354 fprintf(fd,"row %d:",i); 355 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 356 #if defined(PETSC_COMPLEX) 357 if (imag(a->a[j]) != 0.0 && real(a->a[j]) != 0.0) 358 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 359 else if (real(a->a[j]) != 0.0) 360 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 361 #else 362 if (a->a[j] != 0.0) fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 363 #endif 364 } 365 fprintf(fd,"\n"); 366 } 367 } 368 else { 369 for ( i=0; i<m; i++ ) { 370 fprintf(fd,"row %d:",i); 371 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 372 #if defined(PETSC_COMPLEX) 373 if (imag(a->a[j]) != 0.0) { 374 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 375 } 376 else { 377 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 378 } 379 #else 380 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 381 #endif 382 } 383 fprintf(fd,"\n"); 384 } 385 } 386 fflush(fd); 387 return 0; 388 } 389 390 #undef __FUNCTION__ 391 #define __FUNCTION__ "MatView_SeqAIJ_Draw" 392 static int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 393 { 394 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 395 int ierr, i,j, m = a->m, shift = a->indexshift,pause,color; 396 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r; 397 Draw draw; 398 DrawButton button; 399 PetscTruth isnull; 400 401 ViewerDrawGetDraw(viewer,&draw); 402 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) return 0; 403 404 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 405 xr += w; yr += h; xl = -w; yl = -h; 406 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 407 /* loop over matrix elements drawing boxes */ 408 color = DRAW_BLUE; 409 for ( i=0; i<m; i++ ) { 410 y_l = m - i - 1.0; y_r = y_l + 1.0; 411 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 412 x_l = a->j[j] + shift; x_r = x_l + 1.0; 413 #if defined(PETSC_COMPLEX) 414 if (real(a->a[j]) >= 0.) continue; 415 #else 416 if (a->a[j] >= 0.) continue; 417 #endif 418 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 419 } 420 } 421 color = DRAW_CYAN; 422 for ( i=0; i<m; i++ ) { 423 y_l = m - i - 1.0; y_r = y_l + 1.0; 424 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 425 x_l = a->j[j] + shift; x_r = x_l + 1.0; 426 if (a->a[j] != 0.) continue; 427 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 428 } 429 } 430 color = DRAW_RED; 431 for ( i=0; i<m; i++ ) { 432 y_l = m - i - 1.0; y_r = y_l + 1.0; 433 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 434 x_l = a->j[j] + shift; x_r = x_l + 1.0; 435 #if defined(PETSC_COMPLEX) 436 if (real(a->a[j]) <= 0.) continue; 437 #else 438 if (a->a[j] <= 0.) continue; 439 #endif 440 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 441 } 442 } 443 DrawFlush(draw); 444 DrawGetPause(draw,&pause); 445 if (pause >= 0) { PetscSleep(pause); return 0;} 446 447 /* allow the matrix to zoom or shrink */ 448 ierr = DrawCheckResizedWindow(draw); 449 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 450 while (button != BUTTON_RIGHT) { 451 DrawClear(draw); 452 if (button == BUTTON_LEFT) scale = .5; 453 else if (button == BUTTON_CENTER) scale = 2.; 454 xl = scale*(xl + w - xc) + xc - w*scale; 455 xr = scale*(xr - w - xc) + xc + w*scale; 456 yl = scale*(yl + h - yc) + yc - h*scale; 457 yr = scale*(yr - h - yc) + yc + h*scale; 458 w *= scale; h *= scale; 459 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 460 color = DRAW_BLUE; 461 for ( i=0; i<m; i++ ) { 462 y_l = m - i - 1.0; y_r = y_l + 1.0; 463 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 464 x_l = a->j[j] + shift; x_r = x_l + 1.0; 465 #if defined(PETSC_COMPLEX) 466 if (real(a->a[j]) >= 0.) continue; 467 #else 468 if (a->a[j] >= 0.) continue; 469 #endif 470 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 471 } 472 } 473 color = DRAW_CYAN; 474 for ( i=0; i<m; i++ ) { 475 y_l = m - i - 1.0; y_r = y_l + 1.0; 476 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 477 x_l = a->j[j] + shift; x_r = x_l + 1.0; 478 if (a->a[j] != 0.) continue; 479 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 480 } 481 } 482 color = DRAW_RED; 483 for ( i=0; i<m; i++ ) { 484 y_l = m - i - 1.0; y_r = y_l + 1.0; 485 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 486 x_l = a->j[j] + shift; x_r = x_l + 1.0; 487 #if defined(PETSC_COMPLEX) 488 if (real(a->a[j]) <= 0.) continue; 489 #else 490 if (a->a[j] <= 0.) continue; 491 #endif 492 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 493 } 494 } 495 ierr = DrawCheckResizedWindow(draw); 496 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 497 } 498 return 0; 499 } 500 501 #undef __FUNCTION__ 502 #define __FUNCTION__ "MatView_SeqAIJ" 503 static int MatView_SeqAIJ(PetscObject obj,Viewer viewer) 504 { 505 Mat A = (Mat) obj; 506 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 507 ViewerType vtype; 508 int ierr; 509 510 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 511 if (vtype == MATLAB_VIEWER) { 512 return ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); 513 } 514 else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 515 return MatView_SeqAIJ_ASCII(A,viewer); 516 } 517 else if (vtype == BINARY_FILE_VIEWER) { 518 return MatView_SeqAIJ_Binary(A,viewer); 519 } 520 else if (vtype == DRAW_VIEWER) { 521 return MatView_SeqAIJ_Draw(A,viewer); 522 } 523 return 0; 524 } 525 526 extern int Mat_AIJ_CheckInode(Mat); 527 #undef __FUNCTION__ 528 #define __FUNCTION__ "MatAssemblyEnd_SeqAIJ" 529 static int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 530 { 531 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 532 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 533 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift; 534 Scalar *aa = a->a, *ap; 535 536 if (mode == MAT_FLUSH_ASSEMBLY) return 0; 537 538 for ( i=1; i<m; i++ ) { 539 /* move each row back by the amount of empty slots (fshift) before it*/ 540 fshift += imax[i-1] - ailen[i-1]; 541 if (fshift) { 542 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 543 N = ailen[i]; 544 for ( j=0; j<N; j++ ) { 545 ip[j-fshift] = ip[j]; 546 ap[j-fshift] = ap[j]; 547 } 548 } 549 ai[i] = ai[i-1] + ailen[i-1]; 550 } 551 if (m) { 552 fshift += imax[m-1] - ailen[m-1]; 553 ai[m] = ai[m-1] + ailen[m-1]; 554 } 555 /* reset ilen and imax for each row */ 556 for ( i=0; i<m; i++ ) { 557 ailen[i] = imax[i] = ai[i+1] - ai[i]; 558 } 559 a->nz = ai[m] + shift; 560 561 /* diagonals may have moved, so kill the diagonal pointers */ 562 if (fshift && a->diag) { 563 PetscFree(a->diag); 564 PLogObjectMemory(A,-(m+1)*sizeof(int)); 565 a->diag = 0; 566 } 567 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 568 m,a->n,fshift,a->nz); 569 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 570 a->reallocs); 571 a->reallocs = 0; 572 A->info.nz_unneeded = (double)fshift; 573 574 /* check out for identical nodes. If found, use inode functions */ 575 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 576 return 0; 577 } 578 579 #undef __FUNCTION__ 580 #define __FUNCTION__ "MatZeroEntries_SeqAIJ" 581 static int MatZeroEntries_SeqAIJ(Mat A) 582 { 583 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 584 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 585 return 0; 586 } 587 588 #undef __FUNCTION__ 589 #define __FUNCTION__ "MatDestroy_SeqAIJ" 590 int MatDestroy_SeqAIJ(PetscObject obj) 591 { 592 Mat A = (Mat) obj; 593 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 594 int ierr; 595 596 #if defined(PETSC_LOG) 597 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 598 #endif 599 PetscFree(a->a); 600 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 601 if (a->diag) PetscFree(a->diag); 602 if (a->ilen) PetscFree(a->ilen); 603 if (a->imax) PetscFree(a->imax); 604 if (a->solve_work) PetscFree(a->solve_work); 605 if (a->inode.size) PetscFree(a->inode.size); 606 PetscFree(a); 607 if (A->mapping) { 608 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 609 } 610 PLogObjectDestroy(A); 611 PetscHeaderDestroy(A); 612 return 0; 613 } 614 615 #undef __FUNCTION__ 616 #define __FUNCTION__ "MatCompress_SeqAIJ" 617 static int MatCompress_SeqAIJ(Mat A) 618 { 619 return 0; 620 } 621 622 #undef __FUNCTION__ 623 #define __FUNCTION__ "MatSetOption_SeqAIJ" 624 static int MatSetOption_SeqAIJ(Mat A,MatOption op) 625 { 626 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 627 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 628 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 629 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 630 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 631 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 632 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 633 else if (op == MAT_ROWS_SORTED || 634 op == MAT_ROWS_UNSORTED || 635 op == MAT_SYMMETRIC || 636 op == MAT_STRUCTURALLY_SYMMETRIC || 637 op == MAT_YES_NEW_DIAGONALS || 638 op == MAT_IGNORE_OFF_PROCESSOR_ENTRIES) 639 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 640 else if (op == MAT_NO_NEW_DIAGONALS) 641 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:MAT_NO_NEW_DIAGONALS");} 642 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 643 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 644 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 645 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 646 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 647 else 648 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:unknown option");} 649 return 0; 650 } 651 652 #undef __FUNCTION__ 653 #define __FUNCTION__ "MatGetDiagonal_SeqAIJ" 654 static int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 655 { 656 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 657 int i,j, n,shift = a->indexshift; 658 Scalar *x, zero = 0.0; 659 660 VecSet(&zero,v); 661 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 662 if (n != a->m) SETERRQ(1,"MatGetDiagonal_SeqAIJ:Nonconforming matrix and vector"); 663 for ( i=0; i<a->m; i++ ) { 664 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 665 if (a->j[j]+shift == i) { 666 x[i] = a->a[j]; 667 break; 668 } 669 } 670 } 671 return 0; 672 } 673 674 /* -------------------------------------------------------*/ 675 /* Should check that shapes of vectors and matrices match */ 676 /* -------------------------------------------------------*/ 677 #undef __FUNCTION__ 678 #define __FUNCTION__ "MatMultTrans_SeqAIJ" 679 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 680 { 681 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 682 Scalar *x, *y, *v, alpha; 683 int m = a->m, n, i, *idx, shift = a->indexshift; 684 685 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 686 PetscMemzero(y,a->n*sizeof(Scalar)); 687 y = y + shift; /* shift for Fortran start by 1 indexing */ 688 for ( i=0; i<m; i++ ) { 689 idx = a->j + a->i[i] + shift; 690 v = a->a + a->i[i] + shift; 691 n = a->i[i+1] - a->i[i]; 692 alpha = x[i]; 693 while (n-->0) {y[*idx++] += alpha * *v++;} 694 } 695 PLogFlops(2*a->nz - a->n); 696 return 0; 697 } 698 699 #undef __FUNCTION__ 700 #define __FUNCTION__ "MatMultTransAdd_SeqAIJ" 701 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 702 { 703 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 704 Scalar *x, *y, *v, alpha; 705 int m = a->m, n, i, *idx,shift = a->indexshift; 706 707 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 708 if (zz != yy) VecCopy(zz,yy); 709 y = y + shift; /* shift for Fortran start by 1 indexing */ 710 for ( i=0; i<m; i++ ) { 711 idx = a->j + a->i[i] + shift; 712 v = a->a + a->i[i] + shift; 713 n = a->i[i+1] - a->i[i]; 714 alpha = x[i]; 715 while (n-->0) {y[*idx++] += alpha * *v++;} 716 } 717 PLogFlops(2*a->nz); 718 return 0; 719 } 720 721 #undef __FUNCTION__ 722 #define __FUNCTION__ "MatMult_SeqAIJ" 723 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 724 { 725 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 726 Scalar *x, *y, *v, sum; 727 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 728 729 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 730 x = x + shift; /* shift for Fortran start by 1 indexing */ 731 idx = a->j; 732 v = a->a; 733 ii = a->i; 734 for ( i=0; i<m; i++ ) { 735 n = ii[1] - ii[0]; ii++; 736 sum = 0.0; 737 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 738 /* for ( j=n-1; j>-1; j--) sum += v[j]*x[idx[j]]; */ 739 while (n--) sum += *v++ * x[*idx++]; 740 y[i] = sum; 741 } 742 PLogFlops(2*a->nz - m); 743 return 0; 744 } 745 746 #undef __FUNCTION__ 747 #define __FUNCTION__ "MatMultAdd_SeqAIJ" 748 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 749 { 750 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 751 Scalar *x, *y, *z, *v, sum; 752 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 753 754 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 755 x = x + shift; /* shift for Fortran start by 1 indexing */ 756 idx = a->j; 757 v = a->a; 758 ii = a->i; 759 for ( i=0; i<m; i++ ) { 760 n = ii[1] - ii[0]; ii++; 761 sum = y[i]; 762 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 763 while (n--) sum += *v++ * x[*idx++]; 764 z[i] = sum; 765 } 766 PLogFlops(2*a->nz); 767 return 0; 768 } 769 770 /* 771 Adds diagonal pointers to sparse matrix structure. 772 */ 773 774 #undef __FUNCTION__ 775 #define __FUNCTION__ "MatMarkDiag_SeqAIJ" 776 int MatMarkDiag_SeqAIJ(Mat A) 777 { 778 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 779 int i,j, *diag, m = a->m,shift = a->indexshift; 780 781 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 782 PLogObjectMemory(A,(m+1)*sizeof(int)); 783 for ( i=0; i<a->m; i++ ) { 784 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 785 if (a->j[j]+shift == i) { 786 diag[i] = j - shift; 787 break; 788 } 789 } 790 } 791 a->diag = diag; 792 return 0; 793 } 794 795 #undef __FUNCTION__ 796 #define __FUNCTION__ "MatRelax_SeqAIJ" 797 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 798 double fshift,int its,Vec xx) 799 { 800 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 801 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 802 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 803 804 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 805 if (!a->diag) {if ((ierr = MatMarkDiag_SeqAIJ(A))) return ierr;} 806 diag = a->diag; 807 xs = x + shift; /* shifted by one for index start of a or a->j*/ 808 if (flag == SOR_APPLY_UPPER) { 809 /* apply ( U + D/omega) to the vector */ 810 bs = b + shift; 811 for ( i=0; i<m; i++ ) { 812 d = fshift + a->a[diag[i] + shift]; 813 n = a->i[i+1] - diag[i] - 1; 814 idx = a->j + diag[i] + (!shift); 815 v = a->a + diag[i] + (!shift); 816 sum = b[i]*d/omega; 817 SPARSEDENSEDOT(sum,bs,v,idx,n); 818 x[i] = sum; 819 } 820 return 0; 821 } 822 if (flag == SOR_APPLY_LOWER) { 823 SETERRQ(1,"MatRelax_SeqAIJ:SOR_APPLY_LOWER is not done"); 824 } 825 else if (flag & SOR_EISENSTAT) { 826 /* Let A = L + U + D; where L is lower trianglar, 827 U is upper triangular, E is diagonal; This routine applies 828 829 (L + E)^{-1} A (U + E)^{-1} 830 831 to a vector efficiently using Eisenstat's trick. This is for 832 the case of SSOR preconditioner, so E is D/omega where omega 833 is the relaxation factor. 834 */ 835 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 836 scale = (2.0/omega) - 1.0; 837 838 /* x = (E + U)^{-1} b */ 839 for ( i=m-1; i>=0; i-- ) { 840 d = fshift + a->a[diag[i] + shift]; 841 n = a->i[i+1] - diag[i] - 1; 842 idx = a->j + diag[i] + (!shift); 843 v = a->a + diag[i] + (!shift); 844 sum = b[i]; 845 SPARSEDENSEMDOT(sum,xs,v,idx,n); 846 x[i] = omega*(sum/d); 847 } 848 849 /* t = b - (2*E - D)x */ 850 v = a->a; 851 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 852 853 /* t = (E + L)^{-1}t */ 854 ts = t + shift; /* shifted by one for index start of a or a->j*/ 855 diag = a->diag; 856 for ( i=0; i<m; i++ ) { 857 d = fshift + a->a[diag[i]+shift]; 858 n = diag[i] - a->i[i]; 859 idx = a->j + a->i[i] + shift; 860 v = a->a + a->i[i] + shift; 861 sum = t[i]; 862 SPARSEDENSEMDOT(sum,ts,v,idx,n); 863 t[i] = omega*(sum/d); 864 } 865 866 /* x = x + t */ 867 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 868 PetscFree(t); 869 return 0; 870 } 871 if (flag & SOR_ZERO_INITIAL_GUESS) { 872 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 873 for ( i=0; i<m; i++ ) { 874 d = fshift + a->a[diag[i]+shift]; 875 n = diag[i] - a->i[i]; 876 idx = a->j + a->i[i] + shift; 877 v = a->a + a->i[i] + shift; 878 sum = b[i]; 879 SPARSEDENSEMDOT(sum,xs,v,idx,n); 880 x[i] = omega*(sum/d); 881 } 882 xb = x; 883 } 884 else xb = b; 885 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 886 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 887 for ( i=0; i<m; i++ ) { 888 x[i] *= a->a[diag[i]+shift]; 889 } 890 } 891 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 892 for ( i=m-1; i>=0; i-- ) { 893 d = fshift + a->a[diag[i] + shift]; 894 n = a->i[i+1] - diag[i] - 1; 895 idx = a->j + diag[i] + (!shift); 896 v = a->a + diag[i] + (!shift); 897 sum = xb[i]; 898 SPARSEDENSEMDOT(sum,xs,v,idx,n); 899 x[i] = omega*(sum/d); 900 } 901 } 902 its--; 903 } 904 while (its--) { 905 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 906 for ( i=0; i<m; i++ ) { 907 d = fshift + a->a[diag[i]+shift]; 908 n = a->i[i+1] - a->i[i]; 909 idx = a->j + a->i[i] + shift; 910 v = a->a + a->i[i] + shift; 911 sum = b[i]; 912 SPARSEDENSEMDOT(sum,xs,v,idx,n); 913 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 914 } 915 } 916 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 917 for ( i=m-1; i>=0; i-- ) { 918 d = fshift + a->a[diag[i] + shift]; 919 n = a->i[i+1] - a->i[i]; 920 idx = a->j + a->i[i] + shift; 921 v = a->a + a->i[i] + shift; 922 sum = b[i]; 923 SPARSEDENSEMDOT(sum,xs,v,idx,n); 924 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 925 } 926 } 927 } 928 return 0; 929 } 930 931 #undef __FUNCTION__ 932 #define __FUNCTION__ "MatGetInfo_SeqAIJ" 933 static int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 934 { 935 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 936 937 info->rows_global = (double)a->m; 938 info->columns_global = (double)a->n; 939 info->rows_local = (double)a->m; 940 info->columns_local = (double)a->n; 941 info->block_size = 1.0; 942 info->nz_allocated = (double)a->maxnz; 943 info->nz_used = (double)a->nz; 944 info->nz_unneeded = (double)(a->maxnz - a->nz); 945 /* if (info->nz_unneeded != A->info.nz_unneeded) 946 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 947 info->assemblies = (double)A->num_ass; 948 info->mallocs = (double)a->reallocs; 949 info->memory = A->mem; 950 if (A->factor) { 951 info->fill_ratio_given = A->info.fill_ratio_given; 952 info->fill_ratio_needed = A->info.fill_ratio_needed; 953 info->factor_mallocs = A->info.factor_mallocs; 954 } else { 955 info->fill_ratio_given = 0; 956 info->fill_ratio_needed = 0; 957 info->factor_mallocs = 0; 958 } 959 return 0; 960 } 961 962 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 963 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 964 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 965 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 966 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 967 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 968 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 969 970 #undef __FUNCTION__ 971 #define __FUNCTION__ "MatZeroRows_SeqAIJ" 972 static int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 973 { 974 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 975 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 976 977 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 978 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 979 if (diag) { 980 for ( i=0; i<N; i++ ) { 981 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 982 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 983 a->ilen[rows[i]] = 1; 984 a->a[a->i[rows[i]]+shift] = *diag; 985 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 986 } 987 else { 988 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 989 CHKERRQ(ierr); 990 } 991 } 992 } 993 else { 994 for ( i=0; i<N; i++ ) { 995 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 996 a->ilen[rows[i]] = 0; 997 } 998 } 999 ISRestoreIndices(is,&rows); 1000 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1001 return 0; 1002 } 1003 1004 #undef __FUNCTION__ 1005 #define __FUNCTION__ "MatGetSize_SeqAIJ" 1006 static int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1007 { 1008 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1009 *m = a->m; *n = a->n; 1010 return 0; 1011 } 1012 1013 #undef __FUNCTION__ 1014 #define __FUNCTION__ "MatGetOwnershipRange_SeqAIJ" 1015 static int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1016 { 1017 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1018 *m = 0; *n = a->m; 1019 return 0; 1020 } 1021 1022 #undef __FUNCTION__ 1023 #define __FUNCTION__ "MatGetRow_SeqAIJ" 1024 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1025 { 1026 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1027 int *itmp,i,shift = a->indexshift; 1028 1029 if (row < 0 || row >= a->m) SETERRQ(1,"MatGetRow_SeqAIJ:Row out of range"); 1030 1031 *nz = a->i[row+1] - a->i[row]; 1032 if (v) *v = a->a + a->i[row] + shift; 1033 if (idx) { 1034 itmp = a->j + a->i[row] + shift; 1035 if (*nz && shift) { 1036 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1037 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1038 } else if (*nz) { 1039 *idx = itmp; 1040 } 1041 else *idx = 0; 1042 } 1043 return 0; 1044 } 1045 1046 #undef __FUNCTION__ 1047 #define __FUNCTION__ "MatRestoreRow_SeqAIJ" 1048 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1049 { 1050 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1051 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1052 return 0; 1053 } 1054 1055 #undef __FUNCTION__ 1056 #define __FUNCTION__ "MatNorm_SeqAIJ" 1057 static int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1058 { 1059 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1060 Scalar *v = a->a; 1061 double sum = 0.0; 1062 int i, j,shift = a->indexshift; 1063 1064 if (type == NORM_FROBENIUS) { 1065 for (i=0; i<a->nz; i++ ) { 1066 #if defined(PETSC_COMPLEX) 1067 sum += real(conj(*v)*(*v)); v++; 1068 #else 1069 sum += (*v)*(*v); v++; 1070 #endif 1071 } 1072 *norm = sqrt(sum); 1073 } 1074 else if (type == NORM_1) { 1075 double *tmp; 1076 int *jj = a->j; 1077 tmp = (double *) PetscMalloc( a->n*sizeof(double) ); CHKPTRQ(tmp); 1078 PetscMemzero(tmp,a->n*sizeof(double)); 1079 *norm = 0.0; 1080 for ( j=0; j<a->nz; j++ ) { 1081 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1082 } 1083 for ( j=0; j<a->n; j++ ) { 1084 if (tmp[j] > *norm) *norm = tmp[j]; 1085 } 1086 PetscFree(tmp); 1087 } 1088 else if (type == NORM_INFINITY) { 1089 *norm = 0.0; 1090 for ( j=0; j<a->m; j++ ) { 1091 v = a->a + a->i[j] + shift; 1092 sum = 0.0; 1093 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1094 sum += PetscAbsScalar(*v); v++; 1095 } 1096 if (sum > *norm) *norm = sum; 1097 } 1098 } 1099 else { 1100 SETERRQ(1,"MatNorm_SeqAIJ:No support for two norm yet"); 1101 } 1102 return 0; 1103 } 1104 1105 #undef __FUNCTION__ 1106 #define __FUNCTION__ "MatTranspose_SeqAIJ" 1107 static int MatTranspose_SeqAIJ(Mat A,Mat *B) 1108 { 1109 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1110 Mat C; 1111 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1112 int shift = a->indexshift; 1113 Scalar *array = a->a; 1114 1115 if (B == PETSC_NULL && m != a->n) 1116 SETERRQ(1,"MatTranspose_SeqAIJ:Square matrix only for in-place"); 1117 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1118 PetscMemzero(col,(1+a->n)*sizeof(int)); 1119 if (shift) { 1120 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1121 } 1122 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1123 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1124 PetscFree(col); 1125 for ( i=0; i<m; i++ ) { 1126 len = ai[i+1]-ai[i]; 1127 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1128 array += len; aj += len; 1129 } 1130 if (shift) { 1131 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1132 } 1133 1134 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1135 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1136 1137 if (B != PETSC_NULL) { 1138 *B = C; 1139 } else { 1140 /* This isn't really an in-place transpose */ 1141 PetscFree(a->a); 1142 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1143 if (a->diag) PetscFree(a->diag); 1144 if (a->ilen) PetscFree(a->ilen); 1145 if (a->imax) PetscFree(a->imax); 1146 if (a->solve_work) PetscFree(a->solve_work); 1147 if (a->inode.size) PetscFree(a->inode.size); 1148 PetscFree(a); 1149 PetscMemcpy(A,C,sizeof(struct _Mat)); 1150 PetscHeaderDestroy(C); 1151 } 1152 return 0; 1153 } 1154 1155 #undef __FUNCTION__ 1156 #define __FUNCTION__ "MatDiagonalScale_SeqAIJ" 1157 static int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1158 { 1159 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1160 Scalar *l,*r,x,*v; 1161 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1162 1163 if (ll) { 1164 /* The local size is used so that VecMPI can be passed to this routine 1165 by MatDiagonalScale_MPIAIJ */ 1166 VecGetLocalSize_Fast(ll,m); 1167 if (m != a->m) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Left scaling vector wrong length"); 1168 VecGetArray_Fast(ll,l); 1169 v = a->a; 1170 for ( i=0; i<m; i++ ) { 1171 x = l[i]; 1172 M = a->i[i+1] - a->i[i]; 1173 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1174 } 1175 PLogFlops(nz); 1176 } 1177 if (rr) { 1178 VecGetLocalSize_Fast(rr,n); 1179 if (n != a->n) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Right scaling vector wrong length"); 1180 VecGetArray_Fast(rr,r); 1181 v = a->a; jj = a->j; 1182 for ( i=0; i<nz; i++ ) { 1183 (*v++) *= r[*jj++ + shift]; 1184 } 1185 PLogFlops(nz); 1186 } 1187 return 0; 1188 } 1189 1190 #undef __FUNCTION__ 1191 #define __FUNCTION__ "MatGetSubMatrix_SeqAIJ" 1192 static int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1193 { 1194 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1195 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1196 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1197 register int sum,lensi; 1198 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1199 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1200 Scalar *a_new,*mat_a; 1201 Mat C; 1202 1203 ierr = ISSorted(isrow,(PetscTruth*)&i); 1204 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:ISrow is not sorted"); 1205 ierr = ISSorted(iscol,(PetscTruth*)&i); 1206 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:IScol is not sorted"); 1207 1208 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1209 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1210 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1211 1212 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1213 /* special case of contiguous rows */ 1214 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1215 starts = lens + ncols; 1216 /* loop over new rows determining lens and starting points */ 1217 for (i=0; i<nrows; i++) { 1218 kstart = ai[irow[i]]+shift; 1219 kend = kstart + ailen[irow[i]]; 1220 for ( k=kstart; k<kend; k++ ) { 1221 if (aj[k]+shift >= first) { 1222 starts[i] = k; 1223 break; 1224 } 1225 } 1226 sum = 0; 1227 while (k < kend) { 1228 if (aj[k++]+shift >= first+ncols) break; 1229 sum++; 1230 } 1231 lens[i] = sum; 1232 } 1233 /* create submatrix */ 1234 if (scall == MAT_REUSE_MATRIX) { 1235 int n_cols,n_rows; 1236 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1237 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1238 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1239 C = *B; 1240 } 1241 else { 1242 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1243 } 1244 c = (Mat_SeqAIJ*) C->data; 1245 1246 /* loop over rows inserting into submatrix */ 1247 a_new = c->a; 1248 j_new = c->j; 1249 i_new = c->i; 1250 i_new[0] = -shift; 1251 for (i=0; i<nrows; i++) { 1252 ii = starts[i]; 1253 lensi = lens[i]; 1254 for ( k=0; k<lensi; k++ ) { 1255 *j_new++ = aj[ii+k] - first; 1256 } 1257 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1258 a_new += lensi; 1259 i_new[i+1] = i_new[i] + lensi; 1260 c->ilen[i] = lensi; 1261 } 1262 PetscFree(lens); 1263 } 1264 else { 1265 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1266 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1267 ssmap = smap + shift; 1268 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1269 PetscMemzero(smap,oldcols*sizeof(int)); 1270 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1271 /* determine lens of each row */ 1272 for (i=0; i<nrows; i++) { 1273 kstart = ai[irow[i]]+shift; 1274 kend = kstart + a->ilen[irow[i]]; 1275 lens[i] = 0; 1276 for ( k=kstart; k<kend; k++ ) { 1277 if (ssmap[aj[k]]) { 1278 lens[i]++; 1279 } 1280 } 1281 } 1282 /* Create and fill new matrix */ 1283 if (scall == MAT_REUSE_MATRIX) { 1284 c = (Mat_SeqAIJ *)((*B)->data); 1285 1286 if (c->m != nrows || c->n != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1287 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1288 SETERRQ(1,"MatGetSubmatrices_SeqAIJ:Cannot reuse matrix. wrong no of nonzeros"); 1289 } 1290 PetscMemzero(c->ilen,c->m*sizeof(int)); 1291 C = *B; 1292 } 1293 else { 1294 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1295 } 1296 c = (Mat_SeqAIJ *)(C->data); 1297 for (i=0; i<nrows; i++) { 1298 row = irow[i]; 1299 nznew = 0; 1300 kstart = ai[row]+shift; 1301 kend = kstart + a->ilen[row]; 1302 mat_i = c->i[i]+shift; 1303 mat_j = c->j + mat_i; 1304 mat_a = c->a + mat_i; 1305 mat_ilen = c->ilen + i; 1306 for ( k=kstart; k<kend; k++ ) { 1307 if ((tcol=ssmap[a->j[k]])) { 1308 *mat_j++ = tcol - (!shift); 1309 *mat_a++ = a->a[k]; 1310 (*mat_ilen)++; 1311 1312 } 1313 } 1314 } 1315 /* Free work space */ 1316 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1317 PetscFree(smap); PetscFree(lens); 1318 } 1319 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1320 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1321 1322 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1323 *B = C; 1324 return 0; 1325 } 1326 1327 /* 1328 note: This can only work for identity for row and col. It would 1329 be good to check this and otherwise generate an error. 1330 */ 1331 #undef __FUNCTION__ 1332 #define __FUNCTION__ "MatILUFactor_SeqAIJ" 1333 static int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1334 { 1335 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1336 int ierr; 1337 Mat outA; 1338 1339 if (fill != 0) SETERRQ(1,"MatILUFactor_SeqAIJ:Only fill=0 supported"); 1340 1341 outA = inA; 1342 inA->factor = FACTOR_LU; 1343 a->row = row; 1344 a->col = col; 1345 1346 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1347 1348 if (!a->diag) { 1349 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1350 } 1351 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1352 return 0; 1353 } 1354 1355 #include "pinclude/plapack.h" 1356 #undef __FUNCTION__ 1357 #define __FUNCTION__ "MatScale_SeqAIJ" 1358 static int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1359 { 1360 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1361 int one = 1; 1362 BLscal_( &a->nz, alpha, a->a, &one ); 1363 PLogFlops(a->nz); 1364 return 0; 1365 } 1366 1367 #undef __FUNCTION__ 1368 #define __FUNCTION__ "MatGetSubMatrices_SeqAIJ" 1369 static int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1370 Mat **B) 1371 { 1372 int ierr,i; 1373 1374 if (scall == MAT_INITIAL_MATRIX) { 1375 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1376 } 1377 1378 for ( i=0; i<n; i++ ) { 1379 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1380 } 1381 return 0; 1382 } 1383 1384 #undef __FUNCTION__ 1385 #define __FUNCTION__ "MatGetBlockSize_SeqAIJ" 1386 static int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1387 { 1388 *bs = 1; 1389 return 0; 1390 } 1391 1392 #undef __FUNCTION__ 1393 #define __FUNCTION__ "MatIncreaseOverlap_SeqAIJ" 1394 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1395 { 1396 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1397 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1398 int start, end, *ai, *aj; 1399 char *table; 1400 shift = a->indexshift; 1401 m = a->m; 1402 ai = a->i; 1403 aj = a->j+shift; 1404 1405 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqAIJ: illegal overlap value used"); 1406 1407 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1408 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1409 1410 for ( i=0; i<is_max; i++ ) { 1411 /* Initialise the two local arrays */ 1412 isz = 0; 1413 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1414 1415 /* Extract the indices, assume there can be duplicate entries */ 1416 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1417 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1418 1419 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 1420 for ( j=0; j<n ; ++j){ 1421 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1422 } 1423 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1424 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1425 1426 k = 0; 1427 for ( j=0; j<ov; j++){ /* for each overlap*/ 1428 n = isz; 1429 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1430 row = nidx[k]; 1431 start = ai[row]; 1432 end = ai[row+1]; 1433 for ( l = start; l<end ; l++){ 1434 val = aj[l] + shift; 1435 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1436 } 1437 } 1438 } 1439 ierr = ISCreateGeneral(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1440 } 1441 PetscFree(table); 1442 PetscFree(nidx); 1443 return 0; 1444 } 1445 1446 #undef __FUNCTION__ 1447 #define __FUNCTION__ "MatPrintHelp_SeqAIJ" 1448 int MatPrintHelp_SeqAIJ(Mat A) 1449 { 1450 static int called = 0; 1451 MPI_Comm comm = A->comm; 1452 1453 if (called) return 0; else called = 1; 1454 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1455 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>\n"); 1456 PetscPrintf(comm," -mat_aij_oneindex - internal indices begin at 1 instead of the default 0.\n"); 1457 PetscPrintf(comm," -mat_aij_no_inode - Do not use inodes\n"); 1458 PetscPrintf(comm," -mat_aij_inode_limit <limit> - Set inode limit (max limit=5)\n"); 1459 #if defined(HAVE_ESSL) 1460 PetscPrintf(comm," -mat_aij_essl - Use IBM sparse LU factorization and solve.\n"); 1461 #endif 1462 return 0; 1463 } 1464 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1465 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1466 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1467 1468 /* -------------------------------------------------------------------*/ 1469 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1470 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1471 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1472 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1473 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1474 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1475 MatLUFactor_SeqAIJ,0, 1476 MatRelax_SeqAIJ, 1477 MatTranspose_SeqAIJ, 1478 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1479 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1480 0,MatAssemblyEnd_SeqAIJ, 1481 MatCompress_SeqAIJ, 1482 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1483 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1484 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1485 MatILUFactorSymbolic_SeqAIJ,0, 1486 0,0,MatConvert_SeqAIJ, 1487 MatConvertSameType_SeqAIJ,0,0, 1488 MatILUFactor_SeqAIJ,0,0, 1489 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1490 MatGetValues_SeqAIJ,0, 1491 MatPrintHelp_SeqAIJ, 1492 MatScale_SeqAIJ,0,0, 1493 MatILUDTFactor_SeqAIJ, 1494 MatGetBlockSize_SeqAIJ, 1495 MatGetRowIJ_SeqAIJ, 1496 MatRestoreRowIJ_SeqAIJ, 1497 MatGetColumnIJ_SeqAIJ, 1498 MatRestoreColumnIJ_SeqAIJ, 1499 MatFDColoringCreate_SeqAIJ, 1500 MatColoringPatch_SeqAIJ}; 1501 1502 extern int MatUseSuperLU_SeqAIJ(Mat); 1503 extern int MatUseEssl_SeqAIJ(Mat); 1504 extern int MatUseDXML_SeqAIJ(Mat); 1505 1506 #undef __FUNCTION__ 1507 #define __FUNCTION__ "MatCreateSeqAIJ" 1508 /*@C 1509 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1510 (the default parallel PETSc format). For good matrix assembly performance 1511 the user should preallocate the matrix storage by setting the parameter nz 1512 (or the array nzz). By setting these parameters accurately, performance 1513 during matrix assembly can be increased by more than a factor of 50. 1514 1515 Input Parameters: 1516 . comm - MPI communicator, set to MPI_COMM_SELF 1517 . m - number of rows 1518 . n - number of columns 1519 . nz - number of nonzeros per row (same for all rows) 1520 . nzz - array containing the number of nonzeros in the various rows 1521 (possibly different for each row) or PETSC_NULL 1522 1523 Output Parameter: 1524 . A - the matrix 1525 1526 Notes: 1527 The AIJ format (also called the Yale sparse matrix format or 1528 compressed row storage), is fully compatible with standard Fortran 77 1529 storage. That is, the stored row and column indices can begin at 1530 either one (as in Fortran) or zero. See the users' manual for details. 1531 1532 Specify the preallocated storage with either nz or nnz (not both). 1533 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1534 allocation. For large problems you MUST preallocate memory or you 1535 will get TERRIBLE performance, see the users' manual chapter on matrices. 1536 1537 By default, this format uses inodes (identical nodes) when possible, to 1538 improve numerical efficiency of Matrix vector products and solves. We 1539 search for consecutive rows with the same nonzero structure, thereby 1540 reusing matrix information to achieve increased efficiency. 1541 1542 Options Database Keys: 1543 $ -mat_aij_no_inode - Do not use inodes 1544 $ -mat_aij_inode_limit <limit> - Set inode limit. 1545 $ (max limit=5) 1546 $ -mat_aij_oneindex - Internally use indexing starting at 1 1547 $ rather than 0. Note: When calling MatSetValues(), 1548 $ the user still MUST index entries starting at 0! 1549 1550 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1551 @*/ 1552 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1553 { 1554 Mat B; 1555 Mat_SeqAIJ *b; 1556 int i, len, ierr, flg,size; 1557 1558 MPI_Comm_size(comm,&size); 1559 if (size > 1) SETERRQ(1,"MatCreateSeqAIJ:Comm must be of size 1"); 1560 1561 *A = 0; 1562 PetscHeaderCreate(B,_Mat,MAT_COOKIE,MATSEQAIJ,comm); 1563 PLogObjectCreate(B); 1564 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1565 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1566 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1567 B->destroy = MatDestroy_SeqAIJ; 1568 B->view = MatView_SeqAIJ; 1569 B->factor = 0; 1570 B->lupivotthreshold = 1.0; 1571 B->mapping = 0; 1572 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1573 &flg); CHKERRQ(ierr); 1574 b->ilu_preserve_row_sums = PETSC_FALSE; 1575 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1576 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1577 b->row = 0; 1578 b->col = 0; 1579 b->indexshift = 0; 1580 b->reallocs = 0; 1581 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1582 if (flg) b->indexshift = -1; 1583 1584 b->m = m; B->m = m; B->M = m; 1585 b->n = n; B->n = n; B->N = n; 1586 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1587 if (nnz == PETSC_NULL) { 1588 if (nz == PETSC_DEFAULT) nz = 10; 1589 else if (nz <= 0) nz = 1; 1590 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1591 nz = nz*m; 1592 } 1593 else { 1594 nz = 0; 1595 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1596 } 1597 1598 /* allocate the matrix space */ 1599 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1600 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1601 b->j = (int *) (b->a + nz); 1602 PetscMemzero(b->j,nz*sizeof(int)); 1603 b->i = b->j + nz; 1604 b->singlemalloc = 1; 1605 1606 b->i[0] = -b->indexshift; 1607 for (i=1; i<m+1; i++) { 1608 b->i[i] = b->i[i-1] + b->imax[i-1]; 1609 } 1610 1611 /* b->ilen will count nonzeros in each row so far. */ 1612 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1613 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1614 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1615 1616 b->nz = 0; 1617 b->maxnz = nz; 1618 b->sorted = 0; 1619 b->roworiented = 1; 1620 b->nonew = 0; 1621 b->diag = 0; 1622 b->solve_work = 0; 1623 b->spptr = 0; 1624 b->inode.node_count = 0; 1625 b->inode.size = 0; 1626 b->inode.limit = 5; 1627 b->inode.max_limit = 5; 1628 B->info.nz_unneeded = (double)b->maxnz; 1629 1630 *A = B; 1631 1632 /* SuperLU is not currently supported through PETSc */ 1633 #if defined(HAVE_SUPERLU) 1634 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1635 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1636 #endif 1637 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1638 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1639 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1640 if (flg) { 1641 if (!b->indexshift) SETERRQ(1,"MatCreateSeqAIJ:need -mat_aij_oneindex with -mat_aij_dxml"); 1642 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1643 } 1644 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1645 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1646 return 0; 1647 } 1648 1649 #undef __FUNCTION__ 1650 #define __FUNCTION__ "MatConvertSameType_SeqAIJ" 1651 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1652 { 1653 Mat C; 1654 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1655 int i,len, m = a->m,shift = a->indexshift; 1656 1657 *B = 0; 1658 PetscHeaderCreate(C,_Mat,MAT_COOKIE,MATSEQAIJ,A->comm); 1659 PLogObjectCreate(C); 1660 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1661 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1662 C->destroy = MatDestroy_SeqAIJ; 1663 C->view = MatView_SeqAIJ; 1664 C->factor = A->factor; 1665 c->row = 0; 1666 c->col = 0; 1667 c->indexshift = shift; 1668 C->assembled = PETSC_TRUE; 1669 1670 c->m = C->m = a->m; 1671 c->n = C->n = a->n; 1672 C->M = a->m; 1673 C->N = a->n; 1674 1675 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1676 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1677 for ( i=0; i<m; i++ ) { 1678 c->imax[i] = a->imax[i]; 1679 c->ilen[i] = a->ilen[i]; 1680 } 1681 1682 /* allocate the matrix space */ 1683 c->singlemalloc = 1; 1684 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1685 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1686 c->j = (int *) (c->a + a->i[m] + shift); 1687 c->i = c->j + a->i[m] + shift; 1688 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1689 if (m > 0) { 1690 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1691 if (cpvalues == COPY_VALUES) { 1692 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1693 } 1694 } 1695 1696 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1697 c->sorted = a->sorted; 1698 c->roworiented = a->roworiented; 1699 c->nonew = a->nonew; 1700 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1701 1702 if (a->diag) { 1703 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1704 PLogObjectMemory(C,(m+1)*sizeof(int)); 1705 for ( i=0; i<m; i++ ) { 1706 c->diag[i] = a->diag[i]; 1707 } 1708 } 1709 else c->diag = 0; 1710 c->inode.limit = a->inode.limit; 1711 c->inode.max_limit = a->inode.max_limit; 1712 if (a->inode.size){ 1713 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1714 c->inode.node_count = a->inode.node_count; 1715 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1716 } else { 1717 c->inode.size = 0; 1718 c->inode.node_count = 0; 1719 } 1720 c->nz = a->nz; 1721 c->maxnz = a->maxnz; 1722 c->solve_work = 0; 1723 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1724 1725 *B = C; 1726 return 0; 1727 } 1728 1729 #undef __FUNCTION__ 1730 #define __FUNCTION__ "MatLoad_SeqAIJ" 1731 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1732 { 1733 Mat_SeqAIJ *a; 1734 Mat B; 1735 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1736 MPI_Comm comm; 1737 1738 PetscObjectGetComm((PetscObject) viewer,&comm); 1739 MPI_Comm_size(comm,&size); 1740 if (size > 1) SETERRQ(1,"MatLoad_SeqAIJ:view must have one processor"); 1741 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1742 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1743 if (header[0] != MAT_COOKIE) SETERRQ(1,"MatLoad_SeqAIJ:not matrix object in file"); 1744 M = header[1]; N = header[2]; nz = header[3]; 1745 1746 /* read in row lengths */ 1747 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1748 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1749 1750 /* create our matrix */ 1751 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1752 B = *A; 1753 a = (Mat_SeqAIJ *) B->data; 1754 shift = a->indexshift; 1755 1756 /* read in column indices and adjust for Fortran indexing*/ 1757 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1758 if (shift) { 1759 for ( i=0; i<nz; i++ ) { 1760 a->j[i] += 1; 1761 } 1762 } 1763 1764 /* read in nonzero values */ 1765 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1766 1767 /* set matrix "i" values */ 1768 a->i[0] = -shift; 1769 for ( i=1; i<= M; i++ ) { 1770 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1771 a->ilen[i-1] = rowlengths[i-1]; 1772 } 1773 PetscFree(rowlengths); 1774 1775 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1776 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1777 return 0; 1778 } 1779 1780 #undef __FUNCTION__ 1781 #define __FUNCTION__ "MatEqual_SeqAIJ" 1782 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1783 { 1784 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1785 1786 if (B->type !=MATSEQAIJ)SETERRQ(1,"MatEqual_SeqAIJ:Matrices must be same type"); 1787 1788 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1789 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1790 (a->indexshift != b->indexshift)) { 1791 *flg = PETSC_FALSE; return 0; 1792 } 1793 1794 /* if the a->i are the same */ 1795 if (PetscMemcmp(a->i,b->i, (a->n+1)*sizeof(int))) { 1796 *flg = PETSC_FALSE; return 0; 1797 } 1798 1799 /* if a->j are the same */ 1800 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1801 *flg = PETSC_FALSE; return 0; 1802 } 1803 1804 /* if a->a are the same */ 1805 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1806 *flg = PETSC_FALSE; return 0; 1807 } 1808 *flg = PETSC_TRUE; 1809 return 0; 1810 1811 } 1812