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