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