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