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