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