1 2 3 #ifndef lint 4 static char vcid[] = "$Id: aij.c,v 1.207 1997/02/11 19:46:31 curfman Exp curfman $"; 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 __FUNC__ 21 #define __FUNC__ "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,0,"Not implemented"); 28 } 29 30 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 31 32 #undef __FUNC__ 33 #define __FUNC__ "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 __FUNC__ 67 #define __FUNC__ "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 __FUNC__ 83 #define __FUNC__ "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 __FUNC__ 125 #define __FUNC__ "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 __FUNC__ 140 #define __FUNC__ "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,0,"Negative row"); 153 if (row >= a->m) SETERRQ(1,0,"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,0,"Negative column"); 161 if (in[l] >= a->n) SETERRQ(1,0,"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 __FUNC__ 235 #define __FUNC__ "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,0,"Negative row"); 246 if (row >= a->m) SETERRQ(1,0,"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,0,"Negative column"); 251 if (in[l] >= a->n) SETERRQ(1,0,"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 __FUNC__ 278 #define __FUNC__ "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 __FUNC__ 313 #define __FUNC__ "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 __FUNC__ 391 #define __FUNC__ "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 int format; 397 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r,maxv = 0.0; 398 Draw draw; 399 DrawButton button; 400 PetscTruth isnull; 401 402 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 403 ierr = DrawClear(draw); CHKERRQ(ierr); 404 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 405 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) return 0; 406 407 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 408 xr += w; yr += h; xl = -w; yl = -h; 409 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 410 /* loop over matrix elements drawing boxes */ 411 412 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 413 /* Blue for negative, Cyan for zero and Red for positive */ 414 color = DRAW_BLUE; 415 for ( i=0; i<m; i++ ) { 416 y_l = m - i - 1.0; y_r = y_l + 1.0; 417 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 418 x_l = a->j[j] + shift; x_r = x_l + 1.0; 419 #if defined(PETSC_COMPLEX) 420 if (real(a->a[j]) >= 0.) continue; 421 #else 422 if (a->a[j] >= 0.) continue; 423 #endif 424 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 425 } 426 } 427 color = DRAW_CYAN; 428 for ( i=0; i<m; i++ ) { 429 y_l = m - i - 1.0; y_r = y_l + 1.0; 430 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 431 x_l = a->j[j] + shift; x_r = x_l + 1.0; 432 if (a->a[j] != 0.) continue; 433 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 434 } 435 } 436 color = DRAW_RED; 437 for ( i=0; i<m; i++ ) { 438 y_l = m - i - 1.0; y_r = y_l + 1.0; 439 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 440 x_l = a->j[j] + shift; x_r = x_l + 1.0; 441 #if defined(PETSC_COMPLEX) 442 if (real(a->a[j]) <= 0.) continue; 443 #else 444 if (a->a[j] <= 0.) continue; 445 #endif 446 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 447 } 448 } 449 } else { 450 /* use contour shading to indicate magnitude of values */ 451 /* first determine max of all nonzero values */ 452 int nz = a->nz,count; 453 Draw popup; 454 455 for ( i=0; i<nz; i++ ) { 456 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 457 } 458 ierr = DrawCreatePopUp(draw,&popup); CHKERRQ(ierr); 459 ierr = DrawScalePopup(popup,0.0,maxv); CHKERRQ(ierr); 460 count = 0; 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 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 466 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 467 count++; 468 } 469 } 470 } 471 DrawFlush(draw); 472 DrawGetPause(draw,&pause); 473 if (pause >= 0) { PetscSleep(pause); return 0;} 474 475 /* allow the matrix to zoom or shrink */ 476 ierr = DrawCheckResizedWindow(draw); 477 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 478 while (button != BUTTON_RIGHT) { 479 DrawClear(draw); 480 if (button == BUTTON_LEFT) scale = .5; 481 else if (button == BUTTON_CENTER) scale = 2.; 482 xl = scale*(xl + w - xc) + xc - w*scale; 483 xr = scale*(xr - w - xc) + xc + w*scale; 484 yl = scale*(yl + h - yc) + yc - h*scale; 485 yr = scale*(yr - h - yc) + yc + h*scale; 486 w *= scale; h *= scale; 487 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 488 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 489 /* Blue for negative, Cyan for zero and Red for positive */ 490 color = DRAW_BLUE; 491 for ( i=0; i<m; i++ ) { 492 y_l = m - i - 1.0; y_r = y_l + 1.0; 493 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 494 x_l = a->j[j] + shift; x_r = x_l + 1.0; 495 #if defined(PETSC_COMPLEX) 496 if (real(a->a[j]) >= 0.) continue; 497 #else 498 if (a->a[j] >= 0.) continue; 499 #endif 500 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 501 } 502 } 503 color = DRAW_CYAN; 504 for ( i=0; i<m; i++ ) { 505 y_l = m - i - 1.0; y_r = y_l + 1.0; 506 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 507 x_l = a->j[j] + shift; x_r = x_l + 1.0; 508 if (a->a[j] != 0.) continue; 509 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 510 } 511 } 512 color = DRAW_RED; 513 for ( i=0; i<m; i++ ) { 514 y_l = m - i - 1.0; y_r = y_l + 1.0; 515 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 516 x_l = a->j[j] + shift; x_r = x_l + 1.0; 517 #if defined(PETSC_COMPLEX) 518 if (real(a->a[j]) <= 0.) continue; 519 #else 520 if (a->a[j] <= 0.) continue; 521 #endif 522 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 523 } 524 } 525 } else { 526 /* use contour shading to indicate magnitude of values */ 527 int count = 0; 528 for ( i=0; i<m; i++ ) { 529 y_l = m - i - 1.0; y_r = y_l + 1.0; 530 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 531 x_l = a->j[j] + shift; x_r = x_l + 1.0; 532 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 533 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 534 count++; 535 } 536 } 537 } 538 539 ierr = DrawCheckResizedWindow(draw); 540 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 541 } 542 return 0; 543 } 544 545 #undef __FUNC__ 546 #define __FUNC__ "MatView_SeqAIJ" 547 static int MatView_SeqAIJ(PetscObject obj,Viewer viewer) 548 { 549 Mat A = (Mat) obj; 550 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 551 ViewerType vtype; 552 int ierr; 553 554 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 555 if (vtype == MATLAB_VIEWER) { 556 return ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); 557 } 558 else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 559 return MatView_SeqAIJ_ASCII(A,viewer); 560 } 561 else if (vtype == BINARY_FILE_VIEWER) { 562 return MatView_SeqAIJ_Binary(A,viewer); 563 } 564 else if (vtype == DRAW_VIEWER) { 565 return MatView_SeqAIJ_Draw(A,viewer); 566 } 567 return 0; 568 } 569 570 extern int Mat_AIJ_CheckInode(Mat); 571 #undef __FUNC__ 572 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 573 static int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 574 { 575 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 576 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 577 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax; 578 Scalar *aa = a->a, *ap; 579 580 if (mode == MAT_FLUSH_ASSEMBLY) return 0; 581 582 rmax = ailen[0]; /* determine row with most nonzeros */ 583 for ( i=1; i<m; i++ ) { 584 /* move each row back by the amount of empty slots (fshift) before it*/ 585 fshift += imax[i-1] - ailen[i-1]; 586 rmax = PetscMax(rmax,ailen[i]); 587 if (fshift) { 588 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 589 N = ailen[i]; 590 for ( j=0; j<N; j++ ) { 591 ip[j-fshift] = ip[j]; 592 ap[j-fshift] = ap[j]; 593 } 594 } 595 ai[i] = ai[i-1] + ailen[i-1]; 596 } 597 if (m) { 598 fshift += imax[m-1] - ailen[m-1]; 599 ai[m] = ai[m-1] + ailen[m-1]; 600 } 601 /* reset ilen and imax for each row */ 602 for ( i=0; i<m; i++ ) { 603 ailen[i] = imax[i] = ai[i+1] - ai[i]; 604 } 605 a->nz = ai[m] + shift; 606 607 /* diagonals may have moved, so kill the diagonal pointers */ 608 if (fshift && a->diag) { 609 PetscFree(a->diag); 610 PLogObjectMemory(A,-(m+1)*sizeof(int)); 611 a->diag = 0; 612 } 613 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 614 m,a->n,fshift,a->nz); 615 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 616 a->reallocs); 617 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 618 a->reallocs = 0; 619 A->info.nz_unneeded = (double)fshift; 620 621 /* check out for identical nodes. If found, use inode functions */ 622 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 623 return 0; 624 } 625 626 #undef __FUNC__ 627 #define __FUNC__ "MatZeroEntries_SeqAIJ" 628 static int MatZeroEntries_SeqAIJ(Mat A) 629 { 630 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 631 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 632 return 0; 633 } 634 635 #undef __FUNC__ 636 #define __FUNC__ "MatDestroy_SeqAIJ" 637 int MatDestroy_SeqAIJ(PetscObject obj) 638 { 639 Mat A = (Mat) obj; 640 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 641 int ierr; 642 643 #if defined(PETSC_LOG) 644 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 645 #endif 646 PetscFree(a->a); 647 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 648 if (a->diag) PetscFree(a->diag); 649 if (a->ilen) PetscFree(a->ilen); 650 if (a->imax) PetscFree(a->imax); 651 if (a->solve_work) PetscFree(a->solve_work); 652 if (a->inode.size) PetscFree(a->inode.size); 653 PetscFree(a); 654 if (A->mapping) { 655 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 656 } 657 PLogObjectDestroy(A); 658 PetscHeaderDestroy(A); 659 return 0; 660 } 661 662 #undef __FUNC__ 663 #define __FUNC__ "MatCompress_SeqAIJ" 664 static int MatCompress_SeqAIJ(Mat A) 665 { 666 return 0; 667 } 668 669 #undef __FUNC__ 670 #define __FUNC__ "MatSetOption_SeqAIJ" 671 static int MatSetOption_SeqAIJ(Mat A,MatOption op) 672 { 673 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 674 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 675 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 676 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 677 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 678 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 679 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 680 else if (op == MAT_ROWS_SORTED || 681 op == MAT_ROWS_UNSORTED || 682 op == MAT_SYMMETRIC || 683 op == MAT_STRUCTURALLY_SYMMETRIC || 684 op == MAT_YES_NEW_DIAGONALS || 685 op == MAT_IGNORE_OFF_PROCESSOR_ENTRIES) 686 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 687 else if (op == MAT_NO_NEW_DIAGONALS) 688 {SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS");} 689 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 690 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 691 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 692 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 693 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 694 else 695 {SETERRQ(PETSC_ERR_SUP,0,"unknown option");} 696 return 0; 697 } 698 699 #undef __FUNC__ 700 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 701 static int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 702 { 703 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 704 int i,j, n,shift = a->indexshift; 705 Scalar *x, zero = 0.0; 706 707 VecSet(&zero,v); 708 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 709 if (n != a->m) SETERRQ(1,0,"Nonconforming matrix and vector"); 710 for ( i=0; i<a->m; i++ ) { 711 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 712 if (a->j[j]+shift == i) { 713 x[i] = a->a[j]; 714 break; 715 } 716 } 717 } 718 return 0; 719 } 720 721 /* -------------------------------------------------------*/ 722 /* Should check that shapes of vectors and matrices match */ 723 /* -------------------------------------------------------*/ 724 #undef __FUNC__ 725 #define __FUNC__ "MatMultTrans_SeqAIJ" 726 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 727 { 728 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 729 Scalar *x, *y, *v, alpha; 730 int m = a->m, n, i, *idx, shift = a->indexshift; 731 732 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 733 PetscMemzero(y,a->n*sizeof(Scalar)); 734 y = y + shift; /* shift for Fortran start by 1 indexing */ 735 for ( i=0; i<m; i++ ) { 736 idx = a->j + a->i[i] + shift; 737 v = a->a + a->i[i] + shift; 738 n = a->i[i+1] - a->i[i]; 739 alpha = x[i]; 740 while (n-->0) {y[*idx++] += alpha * *v++;} 741 } 742 PLogFlops(2*a->nz - a->n); 743 return 0; 744 } 745 746 #undef __FUNC__ 747 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 748 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 749 { 750 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 751 Scalar *x, *y, *v, alpha; 752 int m = a->m, n, i, *idx,shift = a->indexshift; 753 754 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 755 if (zz != yy) VecCopy(zz,yy); 756 y = y + shift; /* shift for Fortran start by 1 indexing */ 757 for ( i=0; i<m; i++ ) { 758 idx = a->j + a->i[i] + shift; 759 v = a->a + a->i[i] + shift; 760 n = a->i[i+1] - a->i[i]; 761 alpha = x[i]; 762 while (n-->0) {y[*idx++] += alpha * *v++;} 763 } 764 PLogFlops(2*a->nz); 765 return 0; 766 } 767 768 #undef __FUNC__ 769 #define __FUNC__ "MatMult_SeqAIJ" 770 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 771 { 772 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 773 Scalar *x, *y, *v, sum; 774 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 775 776 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 777 x = x + shift; /* shift for Fortran start by 1 indexing */ 778 idx = a->j; 779 v = a->a; 780 ii = a->i; 781 for ( i=0; i<m; i++ ) { 782 n = ii[1] - ii[0]; ii++; 783 sum = 0.0; 784 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 785 /* for ( j=n-1; j>-1; j--) sum += v[j]*x[idx[j]]; */ 786 while (n--) sum += *v++ * x[*idx++]; 787 y[i] = sum; 788 } 789 PLogFlops(2*a->nz - m); 790 return 0; 791 } 792 793 #undef __FUNC__ 794 #define __FUNC__ "MatMultAdd_SeqAIJ" 795 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 796 { 797 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 798 Scalar *x, *y, *z, *v, sum; 799 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 800 801 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 802 x = x + shift; /* shift for Fortran start by 1 indexing */ 803 idx = a->j; 804 v = a->a; 805 ii = a->i; 806 for ( i=0; i<m; i++ ) { 807 n = ii[1] - ii[0]; ii++; 808 sum = y[i]; 809 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 810 while (n--) sum += *v++ * x[*idx++]; 811 z[i] = sum; 812 } 813 PLogFlops(2*a->nz); 814 return 0; 815 } 816 817 /* 818 Adds diagonal pointers to sparse matrix structure. 819 */ 820 821 #undef __FUNC__ 822 #define __FUNC__ "MatMarkDiag_SeqAIJ" 823 int MatMarkDiag_SeqAIJ(Mat A) 824 { 825 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 826 int i,j, *diag, m = a->m,shift = a->indexshift; 827 828 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 829 PLogObjectMemory(A,(m+1)*sizeof(int)); 830 for ( i=0; i<a->m; i++ ) { 831 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 832 if (a->j[j]+shift == i) { 833 diag[i] = j - shift; 834 break; 835 } 836 } 837 } 838 a->diag = diag; 839 return 0; 840 } 841 842 #undef __FUNC__ 843 #define __FUNC__ "MatRelax_SeqAIJ" 844 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 845 double fshift,int its,Vec xx) 846 { 847 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 848 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 849 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 850 851 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 852 if (!a->diag) {if ((ierr = MatMarkDiag_SeqAIJ(A))) return ierr;} 853 diag = a->diag; 854 xs = x + shift; /* shifted by one for index start of a or a->j*/ 855 if (flag == SOR_APPLY_UPPER) { 856 /* apply ( U + D/omega) to the vector */ 857 bs = b + shift; 858 for ( i=0; i<m; i++ ) { 859 d = fshift + a->a[diag[i] + shift]; 860 n = a->i[i+1] - diag[i] - 1; 861 idx = a->j + diag[i] + (!shift); 862 v = a->a + diag[i] + (!shift); 863 sum = b[i]*d/omega; 864 SPARSEDENSEDOT(sum,bs,v,idx,n); 865 x[i] = sum; 866 } 867 return 0; 868 } 869 if (flag == SOR_APPLY_LOWER) { 870 SETERRQ(1,0,"SOR_APPLY_LOWER is not done"); 871 } 872 else if (flag & SOR_EISENSTAT) { 873 /* Let A = L + U + D; where L is lower trianglar, 874 U is upper triangular, E is diagonal; This routine applies 875 876 (L + E)^{-1} A (U + E)^{-1} 877 878 to a vector efficiently using Eisenstat's trick. This is for 879 the case of SSOR preconditioner, so E is D/omega where omega 880 is the relaxation factor. 881 */ 882 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 883 scale = (2.0/omega) - 1.0; 884 885 /* x = (E + U)^{-1} b */ 886 for ( i=m-1; i>=0; i-- ) { 887 d = fshift + a->a[diag[i] + shift]; 888 n = a->i[i+1] - diag[i] - 1; 889 idx = a->j + diag[i] + (!shift); 890 v = a->a + diag[i] + (!shift); 891 sum = b[i]; 892 SPARSEDENSEMDOT(sum,xs,v,idx,n); 893 x[i] = omega*(sum/d); 894 } 895 896 /* t = b - (2*E - D)x */ 897 v = a->a; 898 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 899 900 /* t = (E + L)^{-1}t */ 901 ts = t + shift; /* shifted by one for index start of a or a->j*/ 902 diag = a->diag; 903 for ( i=0; i<m; i++ ) { 904 d = fshift + a->a[diag[i]+shift]; 905 n = diag[i] - a->i[i]; 906 idx = a->j + a->i[i] + shift; 907 v = a->a + a->i[i] + shift; 908 sum = t[i]; 909 SPARSEDENSEMDOT(sum,ts,v,idx,n); 910 t[i] = omega*(sum/d); 911 } 912 913 /* x = x + t */ 914 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 915 PetscFree(t); 916 return 0; 917 } 918 if (flag & SOR_ZERO_INITIAL_GUESS) { 919 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 920 for ( i=0; i<m; i++ ) { 921 d = fshift + a->a[diag[i]+shift]; 922 n = diag[i] - a->i[i]; 923 idx = a->j + a->i[i] + shift; 924 v = a->a + a->i[i] + shift; 925 sum = b[i]; 926 SPARSEDENSEMDOT(sum,xs,v,idx,n); 927 x[i] = omega*(sum/d); 928 } 929 xb = x; 930 } 931 else xb = b; 932 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 933 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 934 for ( i=0; i<m; i++ ) { 935 x[i] *= a->a[diag[i]+shift]; 936 } 937 } 938 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 939 for ( i=m-1; i>=0; i-- ) { 940 d = fshift + a->a[diag[i] + shift]; 941 n = a->i[i+1] - diag[i] - 1; 942 idx = a->j + diag[i] + (!shift); 943 v = a->a + diag[i] + (!shift); 944 sum = xb[i]; 945 SPARSEDENSEMDOT(sum,xs,v,idx,n); 946 x[i] = omega*(sum/d); 947 } 948 } 949 its--; 950 } 951 while (its--) { 952 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 953 for ( i=0; i<m; i++ ) { 954 d = fshift + a->a[diag[i]+shift]; 955 n = a->i[i+1] - a->i[i]; 956 idx = a->j + a->i[i] + shift; 957 v = a->a + a->i[i] + shift; 958 sum = b[i]; 959 SPARSEDENSEMDOT(sum,xs,v,idx,n); 960 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 961 } 962 } 963 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 964 for ( i=m-1; i>=0; i-- ) { 965 d = fshift + a->a[diag[i] + shift]; 966 n = a->i[i+1] - a->i[i]; 967 idx = a->j + a->i[i] + shift; 968 v = a->a + a->i[i] + shift; 969 sum = b[i]; 970 SPARSEDENSEMDOT(sum,xs,v,idx,n); 971 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 972 } 973 } 974 } 975 return 0; 976 } 977 978 #undef __FUNC__ 979 #define __FUNC__ "MatGetInfo_SeqAIJ" 980 static int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 981 { 982 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 983 984 info->rows_global = (double)a->m; 985 info->columns_global = (double)a->n; 986 info->rows_local = (double)a->m; 987 info->columns_local = (double)a->n; 988 info->block_size = 1.0; 989 info->nz_allocated = (double)a->maxnz; 990 info->nz_used = (double)a->nz; 991 info->nz_unneeded = (double)(a->maxnz - a->nz); 992 /* if (info->nz_unneeded != A->info.nz_unneeded) 993 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 994 info->assemblies = (double)A->num_ass; 995 info->mallocs = (double)a->reallocs; 996 info->memory = A->mem; 997 if (A->factor) { 998 info->fill_ratio_given = A->info.fill_ratio_given; 999 info->fill_ratio_needed = A->info.fill_ratio_needed; 1000 info->factor_mallocs = A->info.factor_mallocs; 1001 } else { 1002 info->fill_ratio_given = 0; 1003 info->fill_ratio_needed = 0; 1004 info->factor_mallocs = 0; 1005 } 1006 return 0; 1007 } 1008 1009 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1010 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1011 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1012 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1013 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1014 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1015 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1016 1017 #undef __FUNC__ 1018 #define __FUNC__ "MatZeroRows_SeqAIJ" 1019 static int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1020 { 1021 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1022 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1023 1024 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1025 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1026 if (diag) { 1027 for ( i=0; i<N; i++ ) { 1028 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1029 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1030 a->ilen[rows[i]] = 1; 1031 a->a[a->i[rows[i]]+shift] = *diag; 1032 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1033 } 1034 else { 1035 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 1036 CHKERRQ(ierr); 1037 } 1038 } 1039 } 1040 else { 1041 for ( i=0; i<N; i++ ) { 1042 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1043 a->ilen[rows[i]] = 0; 1044 } 1045 } 1046 ISRestoreIndices(is,&rows); 1047 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1048 return 0; 1049 } 1050 1051 #undef __FUNC__ 1052 #define __FUNC__ "MatGetSize_SeqAIJ" 1053 static int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1054 { 1055 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1056 *m = a->m; *n = a->n; 1057 return 0; 1058 } 1059 1060 #undef __FUNC__ 1061 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1062 static int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1063 { 1064 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1065 *m = 0; *n = a->m; 1066 return 0; 1067 } 1068 1069 #undef __FUNC__ 1070 #define __FUNC__ "MatGetRow_SeqAIJ" 1071 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1072 { 1073 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1074 int *itmp,i,shift = a->indexshift; 1075 1076 if (row < 0 || row >= a->m) SETERRQ(1,0,"Row out of range"); 1077 1078 *nz = a->i[row+1] - a->i[row]; 1079 if (v) *v = a->a + a->i[row] + shift; 1080 if (idx) { 1081 itmp = a->j + a->i[row] + shift; 1082 if (*nz && shift) { 1083 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1084 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1085 } else if (*nz) { 1086 *idx = itmp; 1087 } 1088 else *idx = 0; 1089 } 1090 return 0; 1091 } 1092 1093 #undef __FUNC__ 1094 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1095 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1096 { 1097 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1098 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1099 return 0; 1100 } 1101 1102 #undef __FUNC__ 1103 #define __FUNC__ "MatNorm_SeqAIJ" 1104 static int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1105 { 1106 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1107 Scalar *v = a->a; 1108 double sum = 0.0; 1109 int i, j,shift = a->indexshift; 1110 1111 if (type == NORM_FROBENIUS) { 1112 for (i=0; i<a->nz; i++ ) { 1113 #if defined(PETSC_COMPLEX) 1114 sum += real(conj(*v)*(*v)); v++; 1115 #else 1116 sum += (*v)*(*v); v++; 1117 #endif 1118 } 1119 *norm = sqrt(sum); 1120 } 1121 else if (type == NORM_1) { 1122 double *tmp; 1123 int *jj = a->j; 1124 tmp = (double *) PetscMalloc( a->n*sizeof(double) ); CHKPTRQ(tmp); 1125 PetscMemzero(tmp,a->n*sizeof(double)); 1126 *norm = 0.0; 1127 for ( j=0; j<a->nz; j++ ) { 1128 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1129 } 1130 for ( j=0; j<a->n; j++ ) { 1131 if (tmp[j] > *norm) *norm = tmp[j]; 1132 } 1133 PetscFree(tmp); 1134 } 1135 else if (type == NORM_INFINITY) { 1136 *norm = 0.0; 1137 for ( j=0; j<a->m; j++ ) { 1138 v = a->a + a->i[j] + shift; 1139 sum = 0.0; 1140 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1141 sum += PetscAbsScalar(*v); v++; 1142 } 1143 if (sum > *norm) *norm = sum; 1144 } 1145 } 1146 else { 1147 SETERRQ(1,0,"No support for two norm yet"); 1148 } 1149 return 0; 1150 } 1151 1152 #undef __FUNC__ 1153 #define __FUNC__ "MatTranspose_SeqAIJ" 1154 static int MatTranspose_SeqAIJ(Mat A,Mat *B) 1155 { 1156 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1157 Mat C; 1158 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1159 int shift = a->indexshift; 1160 Scalar *array = a->a; 1161 1162 if (B == PETSC_NULL && m != a->n) 1163 SETERRQ(1,0,"Square matrix only for in-place"); 1164 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1165 PetscMemzero(col,(1+a->n)*sizeof(int)); 1166 if (shift) { 1167 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1168 } 1169 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1170 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1171 PetscFree(col); 1172 for ( i=0; i<m; i++ ) { 1173 len = ai[i+1]-ai[i]; 1174 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1175 array += len; aj += len; 1176 } 1177 if (shift) { 1178 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1179 } 1180 1181 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1182 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1183 1184 if (B != PETSC_NULL) { 1185 *B = C; 1186 } else { 1187 /* This isn't really an in-place transpose */ 1188 PetscFree(a->a); 1189 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1190 if (a->diag) PetscFree(a->diag); 1191 if (a->ilen) PetscFree(a->ilen); 1192 if (a->imax) PetscFree(a->imax); 1193 if (a->solve_work) PetscFree(a->solve_work); 1194 if (a->inode.size) PetscFree(a->inode.size); 1195 PetscFree(a); 1196 PetscMemcpy(A,C,sizeof(struct _Mat)); 1197 PetscHeaderDestroy(C); 1198 } 1199 return 0; 1200 } 1201 1202 #undef __FUNC__ 1203 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1204 static int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1205 { 1206 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1207 Scalar *l,*r,x,*v; 1208 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1209 1210 if (ll) { 1211 /* The local size is used so that VecMPI can be passed to this routine 1212 by MatDiagonalScale_MPIAIJ */ 1213 VecGetLocalSize_Fast(ll,m); 1214 if (m != a->m) SETERRQ(1,0,"Left scaling vector wrong length"); 1215 VecGetArray_Fast(ll,l); 1216 v = a->a; 1217 for ( i=0; i<m; i++ ) { 1218 x = l[i]; 1219 M = a->i[i+1] - a->i[i]; 1220 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1221 } 1222 PLogFlops(nz); 1223 } 1224 if (rr) { 1225 VecGetLocalSize_Fast(rr,n); 1226 if (n != a->n) SETERRQ(1,0,"Right scaling vector wrong length"); 1227 VecGetArray_Fast(rr,r); 1228 v = a->a; jj = a->j; 1229 for ( i=0; i<nz; i++ ) { 1230 (*v++) *= r[*jj++ + shift]; 1231 } 1232 PLogFlops(nz); 1233 } 1234 return 0; 1235 } 1236 1237 #undef __FUNC__ 1238 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1239 static int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1240 { 1241 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1242 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1243 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1244 register int sum,lensi; 1245 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1246 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1247 Scalar *a_new,*mat_a; 1248 Mat C; 1249 1250 ierr = ISSorted(isrow,(PetscTruth*)&i); 1251 if (!i) SETERRQ(1,0,"ISrow is not sorted"); 1252 ierr = ISSorted(iscol,(PetscTruth*)&i); 1253 if (!i) SETERRQ(1,0,"IScol is not sorted"); 1254 1255 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1256 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1257 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1258 1259 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1260 /* special case of contiguous rows */ 1261 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1262 starts = lens + ncols; 1263 /* loop over new rows determining lens and starting points */ 1264 for (i=0; i<nrows; i++) { 1265 kstart = ai[irow[i]]+shift; 1266 kend = kstart + ailen[irow[i]]; 1267 for ( k=kstart; k<kend; k++ ) { 1268 if (aj[k]+shift >= first) { 1269 starts[i] = k; 1270 break; 1271 } 1272 } 1273 sum = 0; 1274 while (k < kend) { 1275 if (aj[k++]+shift >= first+ncols) break; 1276 sum++; 1277 } 1278 lens[i] = sum; 1279 } 1280 /* create submatrix */ 1281 if (scall == MAT_REUSE_MATRIX) { 1282 int n_cols,n_rows; 1283 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1284 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,0,""); 1285 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1286 C = *B; 1287 } 1288 else { 1289 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1290 } 1291 c = (Mat_SeqAIJ*) C->data; 1292 1293 /* loop over rows inserting into submatrix */ 1294 a_new = c->a; 1295 j_new = c->j; 1296 i_new = c->i; 1297 i_new[0] = -shift; 1298 for (i=0; i<nrows; i++) { 1299 ii = starts[i]; 1300 lensi = lens[i]; 1301 for ( k=0; k<lensi; k++ ) { 1302 *j_new++ = aj[ii+k] - first; 1303 } 1304 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1305 a_new += lensi; 1306 i_new[i+1] = i_new[i] + lensi; 1307 c->ilen[i] = lensi; 1308 } 1309 PetscFree(lens); 1310 } 1311 else { 1312 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1313 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1314 ssmap = smap + shift; 1315 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1316 PetscMemzero(smap,oldcols*sizeof(int)); 1317 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1318 /* determine lens of each row */ 1319 for (i=0; i<nrows; i++) { 1320 kstart = ai[irow[i]]+shift; 1321 kend = kstart + a->ilen[irow[i]]; 1322 lens[i] = 0; 1323 for ( k=kstart; k<kend; k++ ) { 1324 if (ssmap[aj[k]]) { 1325 lens[i]++; 1326 } 1327 } 1328 } 1329 /* Create and fill new matrix */ 1330 if (scall == MAT_REUSE_MATRIX) { 1331 c = (Mat_SeqAIJ *)((*B)->data); 1332 1333 if (c->m != nrows || c->n != ncols) SETERRQ(1,0,""); 1334 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1335 SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros"); 1336 } 1337 PetscMemzero(c->ilen,c->m*sizeof(int)); 1338 C = *B; 1339 } 1340 else { 1341 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1342 } 1343 c = (Mat_SeqAIJ *)(C->data); 1344 for (i=0; i<nrows; i++) { 1345 row = irow[i]; 1346 nznew = 0; 1347 kstart = ai[row]+shift; 1348 kend = kstart + a->ilen[row]; 1349 mat_i = c->i[i]+shift; 1350 mat_j = c->j + mat_i; 1351 mat_a = c->a + mat_i; 1352 mat_ilen = c->ilen + i; 1353 for ( k=kstart; k<kend; k++ ) { 1354 if ((tcol=ssmap[a->j[k]])) { 1355 *mat_j++ = tcol - (!shift); 1356 *mat_a++ = a->a[k]; 1357 (*mat_ilen)++; 1358 1359 } 1360 } 1361 } 1362 /* Free work space */ 1363 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1364 PetscFree(smap); PetscFree(lens); 1365 } 1366 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1367 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1368 1369 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1370 *B = C; 1371 return 0; 1372 } 1373 1374 /* 1375 note: This can only work for identity for row and col. It would 1376 be good to check this and otherwise generate an error. 1377 */ 1378 #undef __FUNC__ 1379 #define __FUNC__ "MatILUFactor_SeqAIJ" 1380 static int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1381 { 1382 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1383 int ierr; 1384 Mat outA; 1385 1386 if (fill != 0) SETERRQ(1,0,"Only fill=0 supported"); 1387 1388 outA = inA; 1389 inA->factor = FACTOR_LU; 1390 a->row = row; 1391 a->col = col; 1392 1393 if (!a->solve_work) { /* this matrix may have been factored before */ 1394 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1395 } 1396 1397 if (!a->diag) { 1398 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1399 } 1400 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1401 return 0; 1402 } 1403 1404 #include "pinclude/plapack.h" 1405 #undef __FUNC__ 1406 #define __FUNC__ "MatScale_SeqAIJ" 1407 static int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1408 { 1409 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1410 int one = 1; 1411 BLscal_( &a->nz, alpha, a->a, &one ); 1412 PLogFlops(a->nz); 1413 return 0; 1414 } 1415 1416 #undef __FUNC__ 1417 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1418 static int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1419 Mat **B) 1420 { 1421 int ierr,i; 1422 1423 if (scall == MAT_INITIAL_MATRIX) { 1424 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1425 } 1426 1427 for ( i=0; i<n; i++ ) { 1428 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1429 } 1430 return 0; 1431 } 1432 1433 #undef __FUNC__ 1434 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1435 static int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1436 { 1437 *bs = 1; 1438 return 0; 1439 } 1440 1441 #undef __FUNC__ 1442 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1443 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1444 { 1445 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1446 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1447 int start, end, *ai, *aj; 1448 char *table; 1449 shift = a->indexshift; 1450 m = a->m; 1451 ai = a->i; 1452 aj = a->j+shift; 1453 1454 if (ov < 0) SETERRQ(1,0,"illegal overlap value used"); 1455 1456 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1457 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1458 1459 for ( i=0; i<is_max; i++ ) { 1460 /* Initialize the two local arrays */ 1461 isz = 0; 1462 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1463 1464 /* Extract the indices, assume there can be duplicate entries */ 1465 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1466 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1467 1468 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1469 for ( j=0; j<n ; ++j){ 1470 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1471 } 1472 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1473 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1474 1475 k = 0; 1476 for ( j=0; j<ov; j++){ /* for each overlap */ 1477 n = isz; 1478 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1479 row = nidx[k]; 1480 start = ai[row]; 1481 end = ai[row+1]; 1482 for ( l = start; l<end ; l++){ 1483 val = aj[l] + shift; 1484 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1485 } 1486 } 1487 } 1488 ierr = ISCreateGeneral(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1489 } 1490 PetscFree(table); 1491 PetscFree(nidx); 1492 return 0; 1493 } 1494 1495 /* -------------------------------------------------------------- */ 1496 #undef __FUNC__ 1497 #define __FUNC__ "MatPermute_SeqAIJ" 1498 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1499 { 1500 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1501 Scalar *vwork; 1502 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1503 int *row,*col,*cnew,j,*lens; 1504 1505 ierr = ISGetIndices(rowp,&row); CHKERRQ(ierr); 1506 ierr = ISGetIndices(colp,&col); CHKERRQ(ierr); 1507 1508 /* determine lengths of permuted rows */ 1509 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1510 for (i=0; i<m; i++ ) { 1511 lens[row[i]] = a->i[i+1] - a->i[i]; 1512 } 1513 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1514 PetscFree(lens); 1515 1516 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1517 for (i=0; i<m; i++) { 1518 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1519 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1520 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1521 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1522 } 1523 PetscFree(cnew); 1524 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1525 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1526 ierr = ISRestoreIndices(rowp,&row); CHKERRQ(ierr); 1527 ierr = ISRestoreIndices(colp,&col); CHKERRQ(ierr); 1528 return 0; 1529 } 1530 1531 #undef __FUNC__ 1532 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1533 int MatPrintHelp_SeqAIJ(Mat A) 1534 { 1535 static int called = 0; 1536 MPI_Comm comm = A->comm; 1537 1538 if (called) return 0; else called = 1; 1539 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1540 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1541 PetscPrintf(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1542 PetscPrintf(comm," -mat_aij_no_inode: Do not use inodes\n"); 1543 PetscPrintf(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1544 #if defined(HAVE_ESSL) 1545 PetscPrintf(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1546 #endif 1547 return 0; 1548 } 1549 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1550 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1551 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1552 1553 /* -------------------------------------------------------------------*/ 1554 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1555 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1556 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1557 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1558 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1559 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1560 MatLUFactor_SeqAIJ,0, 1561 MatRelax_SeqAIJ, 1562 MatTranspose_SeqAIJ, 1563 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1564 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1565 0,MatAssemblyEnd_SeqAIJ, 1566 MatCompress_SeqAIJ, 1567 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1568 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1569 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1570 MatILUFactorSymbolic_SeqAIJ,0, 1571 0,0, 1572 MatConvertSameType_SeqAIJ,0,0, 1573 MatILUFactor_SeqAIJ,0,0, 1574 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1575 MatGetValues_SeqAIJ,0, 1576 MatPrintHelp_SeqAIJ, 1577 MatScale_SeqAIJ,0,0, 1578 MatILUDTFactor_SeqAIJ, 1579 MatGetBlockSize_SeqAIJ, 1580 MatGetRowIJ_SeqAIJ, 1581 MatRestoreRowIJ_SeqAIJ, 1582 MatGetColumnIJ_SeqAIJ, 1583 MatRestoreColumnIJ_SeqAIJ, 1584 MatFDColoringCreate_SeqAIJ, 1585 MatColoringPatch_SeqAIJ, 1586 0, 1587 MatPermute_SeqAIJ}; 1588 1589 extern int MatUseSuperLU_SeqAIJ(Mat); 1590 extern int MatUseEssl_SeqAIJ(Mat); 1591 extern int MatUseDXML_SeqAIJ(Mat); 1592 1593 #undef __FUNC__ 1594 #define __FUNC__ "MatCreateSeqAIJ" 1595 /*@C 1596 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1597 (the default parallel PETSc format). For good matrix assembly performance 1598 the user should preallocate the matrix storage by setting the parameter nz 1599 (or the array nzz). By setting these parameters accurately, performance 1600 during matrix assembly can be increased by more than a factor of 50. 1601 1602 Input Parameters: 1603 . comm - MPI communicator, set to MPI_COMM_SELF 1604 . m - number of rows 1605 . n - number of columns 1606 . nz - number of nonzeros per row (same for all rows) 1607 . nzz - array containing the number of nonzeros in the various rows 1608 (possibly different for each row) or PETSC_NULL 1609 1610 Output Parameter: 1611 . A - the matrix 1612 1613 Notes: 1614 The AIJ format (also called the Yale sparse matrix format or 1615 compressed row storage), is fully compatible with standard Fortran 77 1616 storage. That is, the stored row and column indices can begin at 1617 either one (as in Fortran) or zero. See the users' manual for details. 1618 1619 Specify the preallocated storage with either nz or nnz (not both). 1620 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1621 allocation. For large problems you MUST preallocate memory or you 1622 will get TERRIBLE performance, see the users' manual chapter on matrices. 1623 1624 By default, this format uses inodes (identical nodes) when possible, to 1625 improve numerical efficiency of Matrix vector products and solves. We 1626 search for consecutive rows with the same nonzero structure, thereby 1627 reusing matrix information to achieve increased efficiency. 1628 1629 Options Database Keys: 1630 $ -mat_aij_no_inode - Do not use inodes 1631 $ -mat_aij_inode_limit <limit> - Set inode limit. 1632 $ (max limit=5) 1633 $ -mat_aij_oneindex - Internally use indexing starting at 1 1634 $ rather than 0. Note: When calling MatSetValues(), 1635 $ the user still MUST index entries starting at 0! 1636 1637 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1638 @*/ 1639 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1640 { 1641 Mat B; 1642 Mat_SeqAIJ *b; 1643 int i, len, ierr, flg,size; 1644 1645 MPI_Comm_size(comm,&size); 1646 if (size > 1) SETERRQ(1,0,"Comm must be of size 1"); 1647 1648 *A = 0; 1649 PetscHeaderCreate(B,_Mat,MAT_COOKIE,MATSEQAIJ,comm); 1650 PLogObjectCreate(B); 1651 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1652 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1653 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1654 B->destroy = MatDestroy_SeqAIJ; 1655 B->view = MatView_SeqAIJ; 1656 B->factor = 0; 1657 B->lupivotthreshold = 1.0; 1658 B->mapping = 0; 1659 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1660 &flg); CHKERRQ(ierr); 1661 b->ilu_preserve_row_sums = PETSC_FALSE; 1662 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1663 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1664 b->row = 0; 1665 b->col = 0; 1666 b->indexshift = 0; 1667 b->reallocs = 0; 1668 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1669 if (flg) b->indexshift = -1; 1670 1671 b->m = m; B->m = m; B->M = m; 1672 b->n = n; B->n = n; B->N = n; 1673 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1674 if (nnz == PETSC_NULL) { 1675 if (nz == PETSC_DEFAULT) nz = 10; 1676 else if (nz <= 0) nz = 1; 1677 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1678 nz = nz*m; 1679 } 1680 else { 1681 nz = 0; 1682 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1683 } 1684 1685 /* allocate the matrix space */ 1686 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1687 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1688 b->j = (int *) (b->a + nz); 1689 PetscMemzero(b->j,nz*sizeof(int)); 1690 b->i = b->j + nz; 1691 b->singlemalloc = 1; 1692 1693 b->i[0] = -b->indexshift; 1694 for (i=1; i<m+1; i++) { 1695 b->i[i] = b->i[i-1] + b->imax[i-1]; 1696 } 1697 1698 /* b->ilen will count nonzeros in each row so far. */ 1699 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1700 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1701 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1702 1703 b->nz = 0; 1704 b->maxnz = nz; 1705 b->sorted = 0; 1706 b->roworiented = 1; 1707 b->nonew = 0; 1708 b->diag = 0; 1709 b->solve_work = 0; 1710 b->spptr = 0; 1711 b->inode.node_count = 0; 1712 b->inode.size = 0; 1713 b->inode.limit = 5; 1714 b->inode.max_limit = 5; 1715 B->info.nz_unneeded = (double)b->maxnz; 1716 1717 *A = B; 1718 1719 /* SuperLU is not currently supported through PETSc */ 1720 #if defined(HAVE_SUPERLU) 1721 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1722 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1723 #endif 1724 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1725 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1726 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1727 if (flg) { 1728 if (!b->indexshift) SETERRQ(1,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1729 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1730 } 1731 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1732 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1733 return 0; 1734 } 1735 1736 #undef __FUNC__ 1737 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1738 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1739 { 1740 Mat C; 1741 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1742 int i,len, m = a->m,shift = a->indexshift; 1743 1744 *B = 0; 1745 PetscHeaderCreate(C,_Mat,MAT_COOKIE,MATSEQAIJ,A->comm); 1746 PLogObjectCreate(C); 1747 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1748 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1749 C->destroy = MatDestroy_SeqAIJ; 1750 C->view = MatView_SeqAIJ; 1751 C->factor = A->factor; 1752 c->row = 0; 1753 c->col = 0; 1754 c->indexshift = shift; 1755 C->assembled = PETSC_TRUE; 1756 1757 c->m = C->m = a->m; 1758 c->n = C->n = a->n; 1759 C->M = a->m; 1760 C->N = a->n; 1761 1762 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1763 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1764 for ( i=0; i<m; i++ ) { 1765 c->imax[i] = a->imax[i]; 1766 c->ilen[i] = a->ilen[i]; 1767 } 1768 1769 /* allocate the matrix space */ 1770 c->singlemalloc = 1; 1771 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1772 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1773 c->j = (int *) (c->a + a->i[m] + shift); 1774 c->i = c->j + a->i[m] + shift; 1775 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1776 if (m > 0) { 1777 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1778 if (cpvalues == COPY_VALUES) { 1779 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1780 } 1781 } 1782 1783 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1784 c->sorted = a->sorted; 1785 c->roworiented = a->roworiented; 1786 c->nonew = a->nonew; 1787 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1788 1789 if (a->diag) { 1790 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1791 PLogObjectMemory(C,(m+1)*sizeof(int)); 1792 for ( i=0; i<m; i++ ) { 1793 c->diag[i] = a->diag[i]; 1794 } 1795 } 1796 else c->diag = 0; 1797 c->inode.limit = a->inode.limit; 1798 c->inode.max_limit = a->inode.max_limit; 1799 if (a->inode.size){ 1800 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1801 c->inode.node_count = a->inode.node_count; 1802 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1803 } else { 1804 c->inode.size = 0; 1805 c->inode.node_count = 0; 1806 } 1807 c->nz = a->nz; 1808 c->maxnz = a->maxnz; 1809 c->solve_work = 0; 1810 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1811 1812 *B = C; 1813 return 0; 1814 } 1815 1816 #undef __FUNC__ 1817 #define __FUNC__ "MatLoad_SeqAIJ" 1818 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1819 { 1820 Mat_SeqAIJ *a; 1821 Mat B; 1822 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1823 MPI_Comm comm; 1824 1825 PetscObjectGetComm((PetscObject) viewer,&comm); 1826 MPI_Comm_size(comm,&size); 1827 if (size > 1) SETERRQ(1,0,"view must have one processor"); 1828 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1829 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1830 if (header[0] != MAT_COOKIE) SETERRQ(1,0,"not matrix object in file"); 1831 M = header[1]; N = header[2]; nz = header[3]; 1832 1833 /* read in row lengths */ 1834 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1835 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1836 1837 /* create our matrix */ 1838 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1839 B = *A; 1840 a = (Mat_SeqAIJ *) B->data; 1841 shift = a->indexshift; 1842 1843 /* read in column indices and adjust for Fortran indexing*/ 1844 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1845 if (shift) { 1846 for ( i=0; i<nz; i++ ) { 1847 a->j[i] += 1; 1848 } 1849 } 1850 1851 /* read in nonzero values */ 1852 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1853 1854 /* set matrix "i" values */ 1855 a->i[0] = -shift; 1856 for ( i=1; i<= M; i++ ) { 1857 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1858 a->ilen[i-1] = rowlengths[i-1]; 1859 } 1860 PetscFree(rowlengths); 1861 1862 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1863 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1864 return 0; 1865 } 1866 1867 #undef __FUNC__ 1868 #define __FUNC__ "MatEqual_SeqAIJ" 1869 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1870 { 1871 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1872 1873 if (B->type !=MATSEQAIJ)SETERRQ(1,0,"Matrices must be same type"); 1874 1875 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1876 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1877 (a->indexshift != b->indexshift)) { 1878 *flg = PETSC_FALSE; return 0; 1879 } 1880 1881 /* if the a->i are the same */ 1882 if (PetscMemcmp(a->i,b->i, (a->n+1)*sizeof(int))) { 1883 *flg = PETSC_FALSE; return 0; 1884 } 1885 1886 /* if a->j are the same */ 1887 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1888 *flg = PETSC_FALSE; return 0; 1889 } 1890 1891 /* if a->a are the same */ 1892 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1893 *flg = PETSC_FALSE; return 0; 1894 } 1895 *flg = PETSC_TRUE; 1896 return 0; 1897 1898 } 1899