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