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