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