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