1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.251 1998/03/12 23:18:23 bsmith Exp balay $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "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 PetscOps *Abops; 1289 struct _MatOps *Aops; 1290 1291 /* This isn't really an in-place transpose */ 1292 PetscFree(a->a); 1293 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1294 if (a->diag) PetscFree(a->diag); 1295 if (a->ilen) PetscFree(a->ilen); 1296 if (a->imax) PetscFree(a->imax); 1297 if (a->solve_work) PetscFree(a->solve_work); 1298 if (a->inode.size) PetscFree(a->inode.size); 1299 PetscFree(a); 1300 1301 /* 1302 This is horrible, horrible code. We need to keep the 1303 A pointers for the bops and ops but copy everything 1304 else from C. 1305 */ 1306 Abops = A->bops; 1307 Aops = A->ops; 1308 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1309 A->bops = Abops; 1310 A->ops = Aops; 1311 1312 PetscHeaderDestroy(C); 1313 } 1314 PetscFunctionReturn(0); 1315 } 1316 1317 #undef __FUNC__ 1318 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1319 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1320 { 1321 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1322 Scalar *l,*r,x,*v; 1323 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1324 1325 PetscFunctionBegin; 1326 if (ll) { 1327 /* The local size is used so that VecMPI can be passed to this routine 1328 by MatDiagonalScale_MPIAIJ */ 1329 VecGetLocalSize_Fast(ll,m); 1330 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1331 VecGetArray_Fast(ll,l); 1332 v = a->a; 1333 for ( i=0; i<m; i++ ) { 1334 x = l[i]; 1335 M = a->i[i+1] - a->i[i]; 1336 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1337 } 1338 PLogFlops(nz); 1339 } 1340 if (rr) { 1341 VecGetLocalSize_Fast(rr,n); 1342 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1343 VecGetArray_Fast(rr,r); 1344 v = a->a; jj = a->j; 1345 for ( i=0; i<nz; i++ ) { 1346 (*v++) *= r[*jj++ + shift]; 1347 } 1348 PLogFlops(nz); 1349 } 1350 PetscFunctionReturn(0); 1351 } 1352 1353 #undef __FUNC__ 1354 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1355 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1356 { 1357 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1358 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1359 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1360 register int sum,lensi; 1361 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1362 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1363 Scalar *a_new,*mat_a; 1364 Mat C; 1365 1366 PetscFunctionBegin; 1367 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1368 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1369 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1370 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1371 1372 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1373 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1374 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1375 1376 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1377 /* special case of contiguous rows */ 1378 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1379 starts = lens + ncols; 1380 /* loop over new rows determining lens and starting points */ 1381 for (i=0; i<nrows; i++) { 1382 kstart = ai[irow[i]]+shift; 1383 kend = kstart + ailen[irow[i]]; 1384 for ( k=kstart; k<kend; k++ ) { 1385 if (aj[k]+shift >= first) { 1386 starts[i] = k; 1387 break; 1388 } 1389 } 1390 sum = 0; 1391 while (k < kend) { 1392 if (aj[k++]+shift >= first+ncols) break; 1393 sum++; 1394 } 1395 lens[i] = sum; 1396 } 1397 /* create submatrix */ 1398 if (scall == MAT_REUSE_MATRIX) { 1399 int n_cols,n_rows; 1400 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1401 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1402 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1403 C = *B; 1404 } else { 1405 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1406 } 1407 c = (Mat_SeqAIJ*) C->data; 1408 1409 /* loop over rows inserting into submatrix */ 1410 a_new = c->a; 1411 j_new = c->j; 1412 i_new = c->i; 1413 i_new[0] = -shift; 1414 for (i=0; i<nrows; i++) { 1415 ii = starts[i]; 1416 lensi = lens[i]; 1417 for ( k=0; k<lensi; k++ ) { 1418 *j_new++ = aj[ii+k] - first; 1419 } 1420 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1421 a_new += lensi; 1422 i_new[i+1] = i_new[i] + lensi; 1423 c->ilen[i] = lensi; 1424 } 1425 PetscFree(lens); 1426 } else { 1427 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1428 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1429 ssmap = smap + shift; 1430 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1431 PetscMemzero(smap,oldcols*sizeof(int)); 1432 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1433 /* determine lens of each row */ 1434 for (i=0; i<nrows; i++) { 1435 kstart = ai[irow[i]]+shift; 1436 kend = kstart + a->ilen[irow[i]]; 1437 lens[i] = 0; 1438 for ( k=kstart; k<kend; k++ ) { 1439 if (ssmap[aj[k]]) { 1440 lens[i]++; 1441 } 1442 } 1443 } 1444 /* Create and fill new matrix */ 1445 if (scall == MAT_REUSE_MATRIX) { 1446 c = (Mat_SeqAIJ *)((*B)->data); 1447 1448 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1449 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1450 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1451 } 1452 PetscMemzero(c->ilen,c->m*sizeof(int)); 1453 C = *B; 1454 } else { 1455 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1456 } 1457 c = (Mat_SeqAIJ *)(C->data); 1458 for (i=0; i<nrows; i++) { 1459 row = irow[i]; 1460 nznew = 0; 1461 kstart = ai[row]+shift; 1462 kend = kstart + a->ilen[row]; 1463 mat_i = c->i[i]+shift; 1464 mat_j = c->j + mat_i; 1465 mat_a = c->a + mat_i; 1466 mat_ilen = c->ilen + i; 1467 for ( k=kstart; k<kend; k++ ) { 1468 if ((tcol=ssmap[a->j[k]])) { 1469 *mat_j++ = tcol - (!shift); 1470 *mat_a++ = a->a[k]; 1471 (*mat_ilen)++; 1472 1473 } 1474 } 1475 } 1476 /* Free work space */ 1477 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1478 PetscFree(smap); PetscFree(lens); 1479 } 1480 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1481 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1482 1483 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1484 *B = C; 1485 PetscFunctionReturn(0); 1486 } 1487 1488 /* 1489 note: This can only work for identity for row and col. It would 1490 be good to check this and otherwise generate an error. 1491 */ 1492 #undef __FUNC__ 1493 #define __FUNC__ "MatILUFactor_SeqAIJ" 1494 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1495 { 1496 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1497 int ierr; 1498 Mat outA; 1499 1500 PetscFunctionBegin; 1501 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1502 1503 outA = inA; 1504 inA->factor = FACTOR_LU; 1505 a->row = row; 1506 a->col = col; 1507 1508 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1509 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1510 1511 if (!a->solve_work) { /* this matrix may have been factored before */ 1512 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1513 } 1514 1515 if (!a->diag) { 1516 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1517 } 1518 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1519 PetscFunctionReturn(0); 1520 } 1521 1522 #include "pinclude/plapack.h" 1523 #undef __FUNC__ 1524 #define __FUNC__ "MatScale_SeqAIJ" 1525 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1526 { 1527 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1528 int one = 1; 1529 1530 PetscFunctionBegin; 1531 BLscal_( &a->nz, alpha, a->a, &one ); 1532 PLogFlops(a->nz); 1533 PetscFunctionReturn(0); 1534 } 1535 1536 #undef __FUNC__ 1537 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1538 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1539 { 1540 int ierr,i; 1541 1542 PetscFunctionBegin; 1543 if (scall == MAT_INITIAL_MATRIX) { 1544 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1545 } 1546 1547 for ( i=0; i<n; i++ ) { 1548 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1549 } 1550 PetscFunctionReturn(0); 1551 } 1552 1553 #undef __FUNC__ 1554 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1555 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1556 { 1557 PetscFunctionBegin; 1558 *bs = 1; 1559 PetscFunctionReturn(0); 1560 } 1561 1562 #undef __FUNC__ 1563 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1564 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1565 { 1566 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1567 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1568 int start, end, *ai, *aj; 1569 BT table; 1570 1571 PetscFunctionBegin; 1572 shift = a->indexshift; 1573 m = a->m; 1574 ai = a->i; 1575 aj = a->j+shift; 1576 1577 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1578 1579 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1580 ierr = BTCreate(m,table); CHKERRQ(ierr); 1581 1582 for ( i=0; i<is_max; i++ ) { 1583 /* Initialize the two local arrays */ 1584 isz = 0; 1585 BTMemzero(m,table); 1586 1587 /* Extract the indices, assume there can be duplicate entries */ 1588 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1589 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1590 1591 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1592 for ( j=0; j<n ; ++j){ 1593 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1594 } 1595 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1596 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1597 1598 k = 0; 1599 for ( j=0; j<ov; j++){ /* for each overlap */ 1600 n = isz; 1601 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1602 row = nidx[k]; 1603 start = ai[row]; 1604 end = ai[row+1]; 1605 for ( l = start; l<end ; l++){ 1606 val = aj[l] + shift; 1607 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1608 } 1609 } 1610 } 1611 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1612 } 1613 BTDestroy(table); 1614 PetscFree(nidx); 1615 PetscFunctionReturn(0); 1616 } 1617 1618 /* -------------------------------------------------------------- */ 1619 #undef __FUNC__ 1620 #define __FUNC__ "MatPermute_SeqAIJ" 1621 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1622 { 1623 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1624 Scalar *vwork; 1625 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1626 int *row,*col,*cnew,j,*lens; 1627 IS icolp,irowp; 1628 1629 PetscFunctionBegin; 1630 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1631 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1632 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1633 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1634 1635 /* determine lengths of permuted rows */ 1636 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1637 for (i=0; i<m; i++ ) { 1638 lens[row[i]] = a->i[i+1] - a->i[i]; 1639 } 1640 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1641 PetscFree(lens); 1642 1643 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1644 for (i=0; i<m; i++) { 1645 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1646 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1647 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1648 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1649 } 1650 PetscFree(cnew); 1651 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1652 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1653 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1654 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1655 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1656 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 #undef __FUNC__ 1661 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1662 int MatPrintHelp_SeqAIJ(Mat A) 1663 { 1664 static int called = 0; 1665 MPI_Comm comm = A->comm; 1666 1667 PetscFunctionBegin; 1668 if (called) {PetscFunctionReturn(0);} else called = 1; 1669 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1670 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1671 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1672 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1673 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1674 #if defined(HAVE_ESSL) 1675 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1676 #endif 1677 PetscFunctionReturn(0); 1678 } 1679 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1680 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1681 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1682 1683 /* -------------------------------------------------------------------*/ 1684 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1685 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1686 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1687 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1688 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1689 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1690 MatLUFactor_SeqAIJ,0, 1691 MatRelax_SeqAIJ, 1692 MatTranspose_SeqAIJ, 1693 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1694 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1695 0,MatAssemblyEnd_SeqAIJ, 1696 MatCompress_SeqAIJ, 1697 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1698 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1699 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1700 MatILUFactorSymbolic_SeqAIJ,0, 1701 0,0, 1702 MatConvertSameType_SeqAIJ,0,0, 1703 MatILUFactor_SeqAIJ,0,0, 1704 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1705 MatGetValues_SeqAIJ,0, 1706 MatPrintHelp_SeqAIJ, 1707 MatScale_SeqAIJ,0,0, 1708 MatILUDTFactor_SeqAIJ, 1709 MatGetBlockSize_SeqAIJ, 1710 MatGetRowIJ_SeqAIJ, 1711 MatRestoreRowIJ_SeqAIJ, 1712 MatGetColumnIJ_SeqAIJ, 1713 MatRestoreColumnIJ_SeqAIJ, 1714 MatFDColoringCreate_SeqAIJ, 1715 MatColoringPatch_SeqAIJ, 1716 0, 1717 MatPermute_SeqAIJ}; 1718 1719 extern int MatUseSuperLU_SeqAIJ(Mat); 1720 extern int MatUseEssl_SeqAIJ(Mat); 1721 extern int MatUseDXML_SeqAIJ(Mat); 1722 1723 #undef __FUNC__ 1724 #define __FUNC__ "MatCreateSeqAIJ" 1725 /*@C 1726 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1727 (the default parallel PETSc format). For good matrix assembly performance 1728 the user should preallocate the matrix storage by setting the parameter nz 1729 (or the array nzz). By setting these parameters accurately, performance 1730 during matrix assembly can be increased by more than a factor of 50. 1731 1732 Input Parameters: 1733 . comm - MPI communicator, set to PETSC_COMM_SELF 1734 . m - number of rows 1735 . n - number of columns 1736 . nz - number of nonzeros per row (same for all rows) 1737 . nzz - array containing the number of nonzeros in the various rows 1738 (possibly different for each row) or PETSC_NULL 1739 1740 Output Parameter: 1741 . A - the matrix 1742 1743 Notes: 1744 The AIJ format (also called the Yale sparse matrix format or 1745 compressed row storage), is fully compatible with standard Fortran 77 1746 storage. That is, the stored row and column indices can begin at 1747 either one (as in Fortran) or zero. See the users' manual for details. 1748 1749 Specify the preallocated storage with either nz or nnz (not both). 1750 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1751 allocation. For large problems you MUST preallocate memory or you 1752 will get TERRIBLE performance, see the users' manual chapter on matrices. 1753 1754 By default, this format uses inodes (identical nodes) when possible, to 1755 improve numerical efficiency of matrix-vector products and solves. We 1756 search for consecutive rows with the same nonzero structure, thereby 1757 reusing matrix information to achieve increased efficiency. 1758 1759 Options Database Keys: 1760 $ -mat_aij_no_inode - Do not use inodes 1761 $ -mat_aij_inode_limit <limit> - Set inode limit. 1762 $ (max limit=5) 1763 $ -mat_aij_oneindex - Internally use indexing starting at 1 1764 $ rather than 0. Note: When calling MatSetValues(), 1765 $ the user still MUST index entries starting at 0! 1766 1767 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1768 @*/ 1769 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1770 { 1771 Mat B; 1772 Mat_SeqAIJ *b; 1773 int i, len, ierr, flg,size; 1774 1775 PetscFunctionBegin; 1776 MPI_Comm_size(comm,&size); 1777 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 1778 1779 *A = 0; 1780 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1781 PLogObjectCreate(B); 1782 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1783 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1784 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1785 B->destroy = MatDestroy_SeqAIJ; 1786 B->view = MatView_SeqAIJ; 1787 B->factor = 0; 1788 B->lupivotthreshold = 1.0; 1789 B->mapping = 0; 1790 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1791 &flg); CHKERRQ(ierr); 1792 b->ilu_preserve_row_sums = PETSC_FALSE; 1793 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1794 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1795 b->row = 0; 1796 b->col = 0; 1797 b->icol = 0; 1798 b->indexshift = 0; 1799 b->reallocs = 0; 1800 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1801 if (flg) b->indexshift = -1; 1802 1803 b->m = m; B->m = m; B->M = m; 1804 b->n = n; B->n = n; B->N = n; 1805 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1806 if (nnz == PETSC_NULL) { 1807 if (nz == PETSC_DEFAULT) nz = 10; 1808 else if (nz <= 0) nz = 1; 1809 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1810 nz = nz*m; 1811 } else { 1812 nz = 0; 1813 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1814 } 1815 1816 /* allocate the matrix space */ 1817 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1818 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1819 b->j = (int *) (b->a + nz); 1820 PetscMemzero(b->j,nz*sizeof(int)); 1821 b->i = b->j + nz; 1822 b->singlemalloc = 1; 1823 1824 b->i[0] = -b->indexshift; 1825 for (i=1; i<m+1; i++) { 1826 b->i[i] = b->i[i-1] + b->imax[i-1]; 1827 } 1828 1829 /* b->ilen will count nonzeros in each row so far. */ 1830 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1831 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1832 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1833 1834 b->nz = 0; 1835 b->maxnz = nz; 1836 b->sorted = 0; 1837 b->roworiented = 1; 1838 b->nonew = 0; 1839 b->diag = 0; 1840 b->solve_work = 0; 1841 b->spptr = 0; 1842 b->inode.node_count = 0; 1843 b->inode.size = 0; 1844 b->inode.limit = 5; 1845 b->inode.max_limit = 5; 1846 B->info.nz_unneeded = (double)b->maxnz; 1847 1848 *A = B; 1849 1850 /* SuperLU is not currently supported through PETSc */ 1851 #if defined(HAVE_SUPERLU) 1852 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1853 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1854 #endif 1855 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1856 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1857 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1858 if (flg) { 1859 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1860 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1861 } 1862 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1863 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1864 PetscFunctionReturn(0); 1865 } 1866 1867 #undef __FUNC__ 1868 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1869 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1870 { 1871 Mat C; 1872 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1873 int i,len, m = a->m,shift = a->indexshift; 1874 1875 PetscFunctionBegin; 1876 *B = 0; 1877 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1878 PLogObjectCreate(C); 1879 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1880 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1881 C->destroy = MatDestroy_SeqAIJ; 1882 C->view = MatView_SeqAIJ; 1883 C->factor = A->factor; 1884 c->row = 0; 1885 c->col = 0; 1886 c->icol = 0; 1887 c->indexshift = shift; 1888 C->assembled = PETSC_TRUE; 1889 1890 c->m = C->m = a->m; 1891 c->n = C->n = a->n; 1892 C->M = a->m; 1893 C->N = a->n; 1894 1895 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1896 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1897 for ( i=0; i<m; i++ ) { 1898 c->imax[i] = a->imax[i]; 1899 c->ilen[i] = a->ilen[i]; 1900 } 1901 1902 /* allocate the matrix space */ 1903 c->singlemalloc = 1; 1904 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1905 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1906 c->j = (int *) (c->a + a->i[m] + shift); 1907 c->i = c->j + a->i[m] + shift; 1908 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1909 if (m > 0) { 1910 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1911 if (cpvalues == COPY_VALUES) { 1912 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1913 } 1914 } 1915 1916 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1917 c->sorted = a->sorted; 1918 c->roworiented = a->roworiented; 1919 c->nonew = a->nonew; 1920 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1921 1922 if (a->diag) { 1923 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1924 PLogObjectMemory(C,(m+1)*sizeof(int)); 1925 for ( i=0; i<m; i++ ) { 1926 c->diag[i] = a->diag[i]; 1927 } 1928 } else c->diag = 0; 1929 c->inode.limit = a->inode.limit; 1930 c->inode.max_limit = a->inode.max_limit; 1931 if (a->inode.size){ 1932 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1933 c->inode.node_count = a->inode.node_count; 1934 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1935 } else { 1936 c->inode.size = 0; 1937 c->inode.node_count = 0; 1938 } 1939 c->nz = a->nz; 1940 c->maxnz = a->maxnz; 1941 c->solve_work = 0; 1942 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1943 1944 *B = C; 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #undef __FUNC__ 1949 #define __FUNC__ "MatLoad_SeqAIJ" 1950 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1951 { 1952 Mat_SeqAIJ *a; 1953 Mat B; 1954 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1955 MPI_Comm comm; 1956 1957 PetscFunctionBegin; 1958 PetscObjectGetComm((PetscObject) viewer,&comm); 1959 MPI_Comm_size(comm,&size); 1960 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 1961 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1962 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1963 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 1964 M = header[1]; N = header[2]; nz = header[3]; 1965 1966 if (nz < 0) { 1967 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 1968 } 1969 1970 /* read in row lengths */ 1971 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1972 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1973 1974 /* create our matrix */ 1975 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1976 B = *A; 1977 a = (Mat_SeqAIJ *) B->data; 1978 shift = a->indexshift; 1979 1980 /* read in column indices and adjust for Fortran indexing*/ 1981 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 1982 if (shift) { 1983 for ( i=0; i<nz; i++ ) { 1984 a->j[i] += 1; 1985 } 1986 } 1987 1988 /* read in nonzero values */ 1989 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 1990 1991 /* set matrix "i" values */ 1992 a->i[0] = -shift; 1993 for ( i=1; i<= M; i++ ) { 1994 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1995 a->ilen[i-1] = rowlengths[i-1]; 1996 } 1997 PetscFree(rowlengths); 1998 1999 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2000 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 #undef __FUNC__ 2005 #define __FUNC__ "MatEqual_SeqAIJ" 2006 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2007 { 2008 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2009 2010 PetscFunctionBegin; 2011 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2012 2013 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2014 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2015 (a->indexshift != b->indexshift)) { 2016 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2017 } 2018 2019 /* if the a->i are the same */ 2020 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2021 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2022 } 2023 2024 /* if a->j are the same */ 2025 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2026 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2027 } 2028 2029 /* if a->a are the same */ 2030 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2031 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2032 } 2033 *flg = PETSC_TRUE; 2034 PetscFunctionReturn(0); 2035 2036 } 2037