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