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