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