1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.283 1998/09/26 21:57:04 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 "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,PetscReal(a->a[j]),PetscImaginary(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 (PetscImaginary(a->a[j]) > 0.0 && PetscReal(a->a[j]) != 0.0) 373 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 374 else if (PetscImaginary(a->a[j]) < 0.0 && PetscReal(a->a[j]) != 0.0) 375 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 376 else if (PetscReal(a->a[j]) != 0.0) 377 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(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 (PetscImaginary(a->a[j]) != 0.0 || PetscReal(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 (PetscImaginary(a->a[j]) != 0.0 || PetscReal(a->a[j]) != 0.0) 424 fprintf(fd," %18.16e %18.16e ",PetscReal(a->a[j]),PetscImaginary(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 (PetscImaginary(a->a[j]) > 0.0) { 438 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 439 } else if (PetscImaginary(a->a[j]) < 0.0) { 440 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 441 } else { 442 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(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 (PetscReal(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 (PetscReal(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 + (int)(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 (--A->refct > 0) PetscFunctionReturn(0); 663 664 if (A->mapping) { 665 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 666 } 667 if (A->bmapping) { 668 ierr = ISLocalToGlobalMappingDestroy(A->bmapping); CHKERRQ(ierr); 669 } 670 if (A->rmap) { 671 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 672 } 673 if (A->cmap) { 674 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 675 } 676 #if defined(USE_PETSC_LOG) 677 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 678 #endif 679 PetscFree(a->a); 680 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 681 if (a->diag) PetscFree(a->diag); 682 if (a->ilen) PetscFree(a->ilen); 683 if (a->imax) PetscFree(a->imax); 684 if (a->solve_work) PetscFree(a->solve_work); 685 if (a->inode.size) PetscFree(a->inode.size); 686 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 687 if (a->saved_values) PetscFree(a->saved_values); 688 PetscFree(a); 689 690 PLogObjectDestroy(A); 691 PetscHeaderDestroy(A); 692 PetscFunctionReturn(0); 693 } 694 695 #undef __FUNC__ 696 #define __FUNC__ "MatCompress_SeqAIJ" 697 int MatCompress_SeqAIJ(Mat A) 698 { 699 PetscFunctionBegin; 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNC__ 704 #define __FUNC__ "MatSetOption_SeqAIJ" 705 int MatSetOption_SeqAIJ(Mat A,MatOption op) 706 { 707 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 708 709 PetscFunctionBegin; 710 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 711 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 712 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 713 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 714 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 715 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 716 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 717 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 718 else if (op == MAT_ROWS_SORTED || 719 op == MAT_ROWS_UNSORTED || 720 op == MAT_SYMMETRIC || 721 op == MAT_STRUCTURALLY_SYMMETRIC || 722 op == MAT_YES_NEW_DIAGONALS || 723 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 724 op == MAT_USE_HASH_TABLE) 725 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 726 else if (op == MAT_NO_NEW_DIAGONALS) { 727 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 728 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 729 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 730 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 731 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 732 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 733 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 734 PetscFunctionReturn(0); 735 } 736 737 #undef __FUNC__ 738 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 739 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 740 { 741 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 742 int i,j, n,shift = a->indexshift,ierr; 743 Scalar *x, zero = 0.0; 744 745 PetscFunctionBegin; 746 ierr = VecSet(&zero,v);CHKERRQ(ierr); 747 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 748 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 749 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 750 for ( i=0; i<a->m; i++ ) { 751 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 752 if (a->j[j]+shift == i) { 753 x[i] = a->a[j]; 754 break; 755 } 756 } 757 } 758 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 759 PetscFunctionReturn(0); 760 } 761 762 /* -------------------------------------------------------*/ 763 /* Should check that shapes of vectors and matrices match */ 764 /* -------------------------------------------------------*/ 765 #undef __FUNC__ 766 #define __FUNC__ "MatMultTrans_SeqAIJ" 767 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 768 { 769 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 770 Scalar *x, *y, *v, alpha; 771 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 772 773 PetscFunctionBegin; 774 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 775 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 776 PetscMemzero(y,a->n*sizeof(Scalar)); 777 y = y + shift; /* shift for Fortran start by 1 indexing */ 778 for ( i=0; i<m; i++ ) { 779 idx = a->j + a->i[i] + shift; 780 v = a->a + a->i[i] + shift; 781 n = a->i[i+1] - a->i[i]; 782 alpha = x[i]; 783 while (n-->0) {y[*idx++] += alpha * *v++;} 784 } 785 PLogFlops(2*a->nz - a->n); 786 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 787 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNC__ 792 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 793 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 794 { 795 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 796 Scalar *x, *y, *v, alpha; 797 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 798 799 PetscFunctionBegin; 800 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 801 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 802 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 803 y = y + shift; /* shift for Fortran start by 1 indexing */ 804 for ( i=0; i<m; i++ ) { 805 idx = a->j + a->i[i] + shift; 806 v = a->a + a->i[i] + shift; 807 n = a->i[i+1] - a->i[i]; 808 alpha = x[i]; 809 while (n-->0) {y[*idx++] += alpha * *v++;} 810 } 811 PLogFlops(2*a->nz); 812 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 813 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 814 PetscFunctionReturn(0); 815 } 816 817 #undef __FUNC__ 818 #define __FUNC__ "MatMult_SeqAIJ" 819 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 820 { 821 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 822 Scalar *x, *y, *v, sum; 823 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 824 #if !defined(USE_FORTRAN_KERNEL_MULTAIJ) 825 int n, i, jrow,j; 826 #endif 827 828 #if defined(HAVE_PRAGMA_DISJOINT) 829 #pragma disjoint(*x,*y,*v) 830 #endif 831 832 PetscFunctionBegin; 833 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 834 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 835 x = x + shift; /* shift for Fortran start by 1 indexing */ 836 idx = a->j; 837 v = a->a; 838 ii = a->i; 839 #if defined(USE_FORTRAN_KERNEL_MULTAIJ) 840 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 841 #else 842 v += shift; /* shift for Fortran start by 1 indexing */ 843 idx += shift; 844 for ( i=0; i<m; i++ ) { 845 jrow = ii[i]; 846 n = ii[i+1] - jrow; 847 sum = 0.0; 848 for ( j=0; j<n; j++) { 849 sum += v[jrow]*x[idx[jrow]]; jrow++; 850 } 851 y[i] = sum; 852 } 853 #endif 854 PLogFlops(2*a->nz - m); 855 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 856 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 #undef __FUNC__ 861 #define __FUNC__ "MatMultAdd_SeqAIJ" 862 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 863 { 864 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 865 Scalar *x, *y, *z, *v, sum; 866 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 867 #if !defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 868 int n,i,jrow,j; 869 #endif 870 871 PetscFunctionBegin; 872 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 873 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 874 if (zz != yy) { 875 ierr = VecGetArray(zz,&z); CHKERRQ(ierr); 876 } else { 877 z = y; 878 } 879 x = x + shift; /* shift for Fortran start by 1 indexing */ 880 idx = a->j; 881 v = a->a; 882 ii = a->i; 883 #if defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 884 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 885 #else 886 v += shift; /* shift for Fortran start by 1 indexing */ 887 idx += shift; 888 for ( i=0; i<m; i++ ) { 889 jrow = ii[i]; 890 n = ii[i+1] - jrow; 891 sum = y[i]; 892 for ( j=0; j<n; j++) { 893 sum += v[jrow]*x[idx[jrow]]; jrow++; 894 } 895 z[i] = sum; 896 } 897 #endif 898 PLogFlops(2*a->nz); 899 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 900 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 901 if (zz != yy) { 902 ierr = VecRestoreArray(zz,&z); CHKERRQ(ierr); 903 } 904 PetscFunctionReturn(0); 905 } 906 907 /* 908 Adds diagonal pointers to sparse matrix structure. 909 */ 910 911 #undef __FUNC__ 912 #define __FUNC__ "MatMarkDiag_SeqAIJ" 913 int MatMarkDiag_SeqAIJ(Mat A) 914 { 915 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 916 int i,j, *diag, m = a->m,shift = a->indexshift; 917 918 PetscFunctionBegin; 919 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 920 PLogObjectMemory(A,(m+1)*sizeof(int)); 921 for ( i=0; i<a->m; i++ ) { 922 diag[i] = a->i[i+1]; 923 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 924 if (a->j[j]+shift == i) { 925 diag[i] = j - shift; 926 break; 927 } 928 } 929 } 930 a->diag = diag; 931 PetscFunctionReturn(0); 932 } 933 934 #undef __FUNC__ 935 #define __FUNC__ "MatRelax_SeqAIJ" 936 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 937 double fshift,int its,Vec xx) 938 { 939 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 940 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 941 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 942 943 PetscFunctionBegin; 944 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 945 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 946 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 947 diag = a->diag; 948 xs = x + shift; /* shifted by one for index start of a or a->j*/ 949 if (flag == SOR_APPLY_UPPER) { 950 /* apply ( U + D/omega) to the vector */ 951 bs = b + shift; 952 for ( i=0; i<m; i++ ) { 953 d = fshift + a->a[diag[i] + shift]; 954 n = a->i[i+1] - diag[i] - 1; 955 PLogFlops(2*n-1); 956 idx = a->j + diag[i] + (!shift); 957 v = a->a + diag[i] + (!shift); 958 sum = b[i]*d/omega; 959 SPARSEDENSEDOT(sum,bs,v,idx,n); 960 x[i] = sum; 961 } 962 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 963 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 if (flag == SOR_APPLY_LOWER) { 967 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 968 } else if (flag & SOR_EISENSTAT) { 969 /* Let A = L + U + D; where L is lower trianglar, 970 U is upper triangular, E is diagonal; This routine applies 971 972 (L + E)^{-1} A (U + E)^{-1} 973 974 to a vector efficiently using Eisenstat's trick. This is for 975 the case of SSOR preconditioner, so E is D/omega where omega 976 is the relaxation factor. 977 */ 978 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 979 scale = (2.0/omega) - 1.0; 980 981 /* x = (E + U)^{-1} b */ 982 for ( i=m-1; i>=0; i-- ) { 983 d = fshift + a->a[diag[i] + shift]; 984 n = a->i[i+1] - diag[i] - 1; 985 PLogFlops(2*n-1); 986 idx = a->j + diag[i] + (!shift); 987 v = a->a + diag[i] + (!shift); 988 sum = b[i]; 989 SPARSEDENSEMDOT(sum,xs,v,idx,n); 990 x[i] = omega*(sum/d); 991 } 992 993 /* t = b - (2*E - D)x */ 994 v = a->a; 995 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 996 PLogFlops(2*m); 997 998 999 /* t = (E + L)^{-1}t */ 1000 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1001 diag = a->diag; 1002 for ( i=0; i<m; i++ ) { 1003 d = fshift + a->a[diag[i]+shift]; 1004 n = diag[i] - a->i[i]; 1005 PLogFlops(2*n-1); 1006 idx = a->j + a->i[i] + shift; 1007 v = a->a + a->i[i] + shift; 1008 sum = t[i]; 1009 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1010 t[i] = omega*(sum/d); 1011 } 1012 1013 /* x = x + t */ 1014 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1015 PetscFree(t); 1016 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1017 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1018 PetscFunctionReturn(0); 1019 } 1020 if (flag & SOR_ZERO_INITIAL_GUESS) { 1021 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1022 for ( i=0; i<m; i++ ) { 1023 d = fshift + a->a[diag[i]+shift]; 1024 n = diag[i] - a->i[i]; 1025 PLogFlops(2*n-1); 1026 idx = a->j + a->i[i] + shift; 1027 v = a->a + a->i[i] + shift; 1028 sum = b[i]; 1029 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1030 x[i] = omega*(sum/d); 1031 } 1032 xb = x; 1033 } else xb = b; 1034 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1035 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1036 for ( i=0; i<m; i++ ) { 1037 x[i] *= a->a[diag[i]+shift]; 1038 } 1039 PLogFlops(m); 1040 } 1041 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1042 for ( i=m-1; i>=0; i-- ) { 1043 d = fshift + a->a[diag[i] + shift]; 1044 n = a->i[i+1] - diag[i] - 1; 1045 PLogFlops(2*n-1); 1046 idx = a->j + diag[i] + (!shift); 1047 v = a->a + diag[i] + (!shift); 1048 sum = xb[i]; 1049 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1050 x[i] = omega*(sum/d); 1051 } 1052 } 1053 its--; 1054 } 1055 while (its--) { 1056 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1057 for ( i=0; i<m; i++ ) { 1058 d = fshift + a->a[diag[i]+shift]; 1059 n = a->i[i+1] - a->i[i]; 1060 PLogFlops(2*n-1); 1061 idx = a->j + a->i[i] + shift; 1062 v = a->a + a->i[i] + shift; 1063 sum = b[i]; 1064 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1065 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1066 } 1067 } 1068 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1069 for ( i=m-1; i>=0; i-- ) { 1070 d = fshift + a->a[diag[i] + shift]; 1071 n = a->i[i+1] - a->i[i]; 1072 PLogFlops(2*n-1); 1073 idx = a->j + a->i[i] + shift; 1074 v = a->a + a->i[i] + shift; 1075 sum = b[i]; 1076 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1077 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1078 } 1079 } 1080 } 1081 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1082 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNC__ 1087 #define __FUNC__ "MatGetInfo_SeqAIJ" 1088 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1089 { 1090 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1091 1092 PetscFunctionBegin; 1093 info->rows_global = (double)a->m; 1094 info->columns_global = (double)a->n; 1095 info->rows_local = (double)a->m; 1096 info->columns_local = (double)a->n; 1097 info->block_size = 1.0; 1098 info->nz_allocated = (double)a->maxnz; 1099 info->nz_used = (double)a->nz; 1100 info->nz_unneeded = (double)(a->maxnz - a->nz); 1101 info->assemblies = (double)A->num_ass; 1102 info->mallocs = (double)a->reallocs; 1103 info->memory = A->mem; 1104 if (A->factor) { 1105 info->fill_ratio_given = A->info.fill_ratio_given; 1106 info->fill_ratio_needed = A->info.fill_ratio_needed; 1107 info->factor_mallocs = A->info.factor_mallocs; 1108 } else { 1109 info->fill_ratio_given = 0; 1110 info->fill_ratio_needed = 0; 1111 info->factor_mallocs = 0; 1112 } 1113 PetscFunctionReturn(0); 1114 } 1115 1116 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1117 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1118 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1119 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1120 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1121 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1122 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1123 1124 #undef __FUNC__ 1125 #define __FUNC__ "MatZeroRows_SeqAIJ" 1126 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1127 { 1128 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1129 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1130 1131 PetscFunctionBegin; 1132 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1133 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1134 if (diag) { 1135 for ( i=0; i<N; i++ ) { 1136 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1137 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1138 a->ilen[rows[i]] = 1; 1139 a->a[a->i[rows[i]]+shift] = *diag; 1140 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1141 } else { 1142 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1143 } 1144 } 1145 } else { 1146 for ( i=0; i<N; i++ ) { 1147 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1148 a->ilen[rows[i]] = 0; 1149 } 1150 } 1151 ISRestoreIndices(is,&rows); 1152 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 #undef __FUNC__ 1157 #define __FUNC__ "MatGetSize_SeqAIJ" 1158 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1159 { 1160 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1161 1162 PetscFunctionBegin; 1163 if (m) *m = a->m; 1164 if (n) *n = a->n; 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNC__ 1169 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1170 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1171 { 1172 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1173 1174 PetscFunctionBegin; 1175 *m = 0; *n = a->m; 1176 PetscFunctionReturn(0); 1177 } 1178 1179 #undef __FUNC__ 1180 #define __FUNC__ "MatGetRow_SeqAIJ" 1181 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1182 { 1183 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1184 int *itmp,i,shift = a->indexshift; 1185 1186 PetscFunctionBegin; 1187 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1188 1189 *nz = a->i[row+1] - a->i[row]; 1190 if (v) *v = a->a + a->i[row] + shift; 1191 if (idx) { 1192 itmp = a->j + a->i[row] + shift; 1193 if (*nz && shift) { 1194 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1195 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1196 } else if (*nz) { 1197 *idx = itmp; 1198 } 1199 else *idx = 0; 1200 } 1201 PetscFunctionReturn(0); 1202 } 1203 1204 #undef __FUNC__ 1205 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1206 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1207 { 1208 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1209 1210 PetscFunctionBegin; 1211 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1212 PetscFunctionReturn(0); 1213 } 1214 1215 #undef __FUNC__ 1216 #define __FUNC__ "MatNorm_SeqAIJ" 1217 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1218 { 1219 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1220 Scalar *v = a->a; 1221 double sum = 0.0; 1222 int i, j,shift = a->indexshift; 1223 1224 PetscFunctionBegin; 1225 if (type == NORM_FROBENIUS) { 1226 for (i=0; i<a->nz; i++ ) { 1227 #if defined(USE_PETSC_COMPLEX) 1228 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1229 #else 1230 sum += (*v)*(*v); v++; 1231 #endif 1232 } 1233 *norm = sqrt(sum); 1234 } else if (type == NORM_1) { 1235 double *tmp; 1236 int *jj = a->j; 1237 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1238 PetscMemzero(tmp,a->n*sizeof(double)); 1239 *norm = 0.0; 1240 for ( j=0; j<a->nz; j++ ) { 1241 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1242 } 1243 for ( j=0; j<a->n; j++ ) { 1244 if (tmp[j] > *norm) *norm = tmp[j]; 1245 } 1246 PetscFree(tmp); 1247 } else if (type == NORM_INFINITY) { 1248 *norm = 0.0; 1249 for ( j=0; j<a->m; j++ ) { 1250 v = a->a + a->i[j] + shift; 1251 sum = 0.0; 1252 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1253 sum += PetscAbsScalar(*v); v++; 1254 } 1255 if (sum > *norm) *norm = sum; 1256 } 1257 } else { 1258 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1259 } 1260 PetscFunctionReturn(0); 1261 } 1262 1263 #undef __FUNC__ 1264 #define __FUNC__ "MatTranspose_SeqAIJ" 1265 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1266 { 1267 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1268 Mat C; 1269 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1270 int shift = a->indexshift; 1271 Scalar *array = a->a; 1272 1273 PetscFunctionBegin; 1274 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1275 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1276 PetscMemzero(col,(1+a->n)*sizeof(int)); 1277 if (shift) { 1278 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1279 } 1280 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1281 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1282 PetscFree(col); 1283 for ( i=0; i<m; i++ ) { 1284 len = ai[i+1]-ai[i]; 1285 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1286 array += len; aj += len; 1287 } 1288 if (shift) { 1289 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1290 } 1291 1292 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1293 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1294 1295 if (B != PETSC_NULL) { 1296 *B = C; 1297 } else { 1298 PetscOps *Abops; 1299 MatOps Aops; 1300 1301 /* This isn't really an in-place transpose */ 1302 PetscFree(a->a); 1303 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1304 if (a->diag) PetscFree(a->diag); 1305 if (a->ilen) PetscFree(a->ilen); 1306 if (a->imax) PetscFree(a->imax); 1307 if (a->solve_work) PetscFree(a->solve_work); 1308 if (a->inode.size) PetscFree(a->inode.size); 1309 PetscFree(a); 1310 1311 1312 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1313 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1314 1315 /* 1316 This is horrible, horrible code. We need to keep the 1317 the bops and ops Structures, copy everything from C 1318 including the function pointers.. 1319 */ 1320 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1321 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1322 Abops = A->bops; 1323 Aops = A->ops; 1324 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1325 A->bops = Abops; 1326 A->ops = Aops; 1327 A->qlist = 0; 1328 1329 PetscHeaderDestroy(C); 1330 } 1331 PetscFunctionReturn(0); 1332 } 1333 1334 #undef __FUNC__ 1335 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1336 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1337 { 1338 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1339 Scalar *l,*r,x,*v; 1340 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1341 1342 PetscFunctionBegin; 1343 if (ll) { 1344 /* The local size is used so that VecMPI can be passed to this routine 1345 by MatDiagonalScale_MPIAIJ */ 1346 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1347 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1348 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1349 v = a->a; 1350 for ( i=0; i<m; i++ ) { 1351 x = l[i]; 1352 M = a->i[i+1] - a->i[i]; 1353 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1354 } 1355 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1356 PLogFlops(nz); 1357 } 1358 if (rr) { 1359 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1360 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1361 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1362 v = a->a; jj = a->j; 1363 for ( i=0; i<nz; i++ ) { 1364 (*v++) *= r[*jj++ + shift]; 1365 } 1366 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1367 PLogFlops(nz); 1368 } 1369 PetscFunctionReturn(0); 1370 } 1371 1372 #undef __FUNC__ 1373 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1374 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1375 { 1376 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1377 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1378 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1379 register int sum,lensi; 1380 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1381 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1382 Scalar *a_new,*mat_a; 1383 Mat C; 1384 PetscTruth stride; 1385 1386 PetscFunctionBegin; 1387 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1388 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1389 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1390 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1391 1392 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1393 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1394 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1395 1396 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1397 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1398 if (stride && step == 1) { 1399 /* special case of contiguous rows */ 1400 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1401 starts = lens + ncols; 1402 /* loop over new rows determining lens and starting points */ 1403 for (i=0; i<nrows; i++) { 1404 kstart = ai[irow[i]]+shift; 1405 kend = kstart + ailen[irow[i]]; 1406 for ( k=kstart; k<kend; k++ ) { 1407 if (aj[k]+shift >= first) { 1408 starts[i] = k; 1409 break; 1410 } 1411 } 1412 sum = 0; 1413 while (k < kend) { 1414 if (aj[k++]+shift >= first+ncols) break; 1415 sum++; 1416 } 1417 lens[i] = sum; 1418 } 1419 /* create submatrix */ 1420 if (scall == MAT_REUSE_MATRIX) { 1421 int n_cols,n_rows; 1422 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1423 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1424 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1425 C = *B; 1426 } else { 1427 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1428 } 1429 c = (Mat_SeqAIJ*) C->data; 1430 1431 /* loop over rows inserting into submatrix */ 1432 a_new = c->a; 1433 j_new = c->j; 1434 i_new = c->i; 1435 i_new[0] = -shift; 1436 for (i=0; i<nrows; i++) { 1437 ii = starts[i]; 1438 lensi = lens[i]; 1439 for ( k=0; k<lensi; k++ ) { 1440 *j_new++ = aj[ii+k] - first; 1441 } 1442 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1443 a_new += lensi; 1444 i_new[i+1] = i_new[i] + lensi; 1445 c->ilen[i] = lensi; 1446 } 1447 PetscFree(lens); 1448 } else { 1449 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1450 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1451 ssmap = smap + shift; 1452 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1453 PetscMemzero(smap,oldcols*sizeof(int)); 1454 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1455 /* determine lens of each row */ 1456 for (i=0; i<nrows; i++) { 1457 kstart = ai[irow[i]]+shift; 1458 kend = kstart + a->ilen[irow[i]]; 1459 lens[i] = 0; 1460 for ( k=kstart; k<kend; k++ ) { 1461 if (ssmap[aj[k]]) { 1462 lens[i]++; 1463 } 1464 } 1465 } 1466 /* Create and fill new matrix */ 1467 if (scall == MAT_REUSE_MATRIX) { 1468 c = (Mat_SeqAIJ *)((*B)->data); 1469 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1470 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1471 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1472 } 1473 PetscMemzero(c->ilen,c->m*sizeof(int)); 1474 C = *B; 1475 } else { 1476 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1477 } 1478 c = (Mat_SeqAIJ *)(C->data); 1479 for (i=0; i<nrows; i++) { 1480 row = irow[i]; 1481 kstart = ai[row]+shift; 1482 kend = kstart + a->ilen[row]; 1483 mat_i = c->i[i]+shift; 1484 mat_j = c->j + mat_i; 1485 mat_a = c->a + mat_i; 1486 mat_ilen = c->ilen + i; 1487 for ( k=kstart; k<kend; k++ ) { 1488 if ((tcol=ssmap[a->j[k]])) { 1489 *mat_j++ = tcol - (!shift); 1490 *mat_a++ = a->a[k]; 1491 (*mat_ilen)++; 1492 1493 } 1494 } 1495 } 1496 /* Free work space */ 1497 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1498 PetscFree(smap); PetscFree(lens); 1499 } 1500 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1501 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1502 1503 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1504 *B = C; 1505 PetscFunctionReturn(0); 1506 } 1507 1508 /* 1509 note: This can only work for identity for row and col. It would 1510 be good to check this and otherwise generate an error. 1511 */ 1512 #undef __FUNC__ 1513 #define __FUNC__ "MatILUFactor_SeqAIJ" 1514 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1515 { 1516 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1517 int ierr; 1518 Mat outA; 1519 1520 PetscFunctionBegin; 1521 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1522 1523 outA = inA; 1524 inA->factor = FACTOR_LU; 1525 a->row = row; 1526 a->col = col; 1527 1528 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1529 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1530 PLogObjectParent(inA,a->icol); 1531 1532 if (!a->solve_work) { /* this matrix may have been factored before */ 1533 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1534 } 1535 1536 if (!a->diag) { 1537 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1538 } 1539 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1540 PetscFunctionReturn(0); 1541 } 1542 1543 #include "pinclude/blaslapack.h" 1544 #undef __FUNC__ 1545 #define __FUNC__ "MatScale_SeqAIJ" 1546 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1547 { 1548 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1549 int one = 1; 1550 1551 PetscFunctionBegin; 1552 BLscal_( &a->nz, alpha, a->a, &one ); 1553 PLogFlops(a->nz); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNC__ 1558 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1559 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1560 { 1561 int ierr,i; 1562 1563 PetscFunctionBegin; 1564 if (scall == MAT_INITIAL_MATRIX) { 1565 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1566 } 1567 1568 for ( i=0; i<n; i++ ) { 1569 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1570 } 1571 PetscFunctionReturn(0); 1572 } 1573 1574 #undef __FUNC__ 1575 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1576 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1577 { 1578 PetscFunctionBegin; 1579 *bs = 1; 1580 PetscFunctionReturn(0); 1581 } 1582 1583 #undef __FUNC__ 1584 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1585 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1586 { 1587 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1588 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1589 int start, end, *ai, *aj; 1590 BT table; 1591 1592 PetscFunctionBegin; 1593 shift = a->indexshift; 1594 m = a->m; 1595 ai = a->i; 1596 aj = a->j+shift; 1597 1598 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1599 1600 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1601 ierr = BTCreate(m,table); CHKERRQ(ierr); 1602 1603 for ( i=0; i<is_max; i++ ) { 1604 /* Initialize the two local arrays */ 1605 isz = 0; 1606 BTMemzero(m,table); 1607 1608 /* Extract the indices, assume there can be duplicate entries */ 1609 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1610 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1611 1612 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1613 for ( j=0; j<n ; ++j){ 1614 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1615 } 1616 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1617 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1618 1619 k = 0; 1620 for ( j=0; j<ov; j++){ /* for each overlap */ 1621 n = isz; 1622 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1623 row = nidx[k]; 1624 start = ai[row]; 1625 end = ai[row+1]; 1626 for ( l = start; l<end ; l++){ 1627 val = aj[l] + shift; 1628 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1629 } 1630 } 1631 } 1632 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1633 } 1634 BTDestroy(table); 1635 PetscFree(nidx); 1636 PetscFunctionReturn(0); 1637 } 1638 1639 /* -------------------------------------------------------------- */ 1640 #undef __FUNC__ 1641 #define __FUNC__ "MatPermute_SeqAIJ" 1642 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1643 { 1644 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1645 Scalar *vwork; 1646 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1647 int *row,*col,*cnew,j,*lens; 1648 IS icolp,irowp; 1649 1650 PetscFunctionBegin; 1651 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1652 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1653 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1654 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1655 1656 /* determine lengths of permuted rows */ 1657 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1658 for (i=0; i<m; i++ ) { 1659 lens[row[i]] = a->i[i+1] - a->i[i]; 1660 } 1661 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1662 PetscFree(lens); 1663 1664 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1665 for (i=0; i<m; i++) { 1666 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1667 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1668 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1669 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1670 } 1671 PetscFree(cnew); 1672 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1673 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1674 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1675 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1676 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1677 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #undef __FUNC__ 1682 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1683 int MatPrintHelp_SeqAIJ(Mat A) 1684 { 1685 static int called = 0; 1686 MPI_Comm comm = A->comm; 1687 1688 PetscFunctionBegin; 1689 if (called) {PetscFunctionReturn(0);} else called = 1; 1690 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1691 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1692 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1693 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1694 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1695 #if defined(HAVE_ESSL) 1696 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1697 #endif 1698 PetscFunctionReturn(0); 1699 } 1700 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1701 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1702 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1703 1704 #undef __FUNC__ 1705 #define __FUNC__ "MatCopy_SeqAIJ" 1706 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1707 { 1708 int ierr; 1709 1710 PetscFunctionBegin; 1711 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1712 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1713 Mat_SeqAIJ *b = (Mat_SeqAIJ *) B->data; 1714 1715 if (a->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1716 if (b->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(B,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1717 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1718 SETERRQ(1,1,"Number of nonzeros in two matrices are different"); 1719 } 1720 PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 1721 } else { 1722 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1723 } 1724 PetscFunctionReturn(0); 1725 } 1726 1727 /* -------------------------------------------------------------------*/ 1728 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1729 MatGetRow_SeqAIJ, 1730 MatRestoreRow_SeqAIJ, 1731 MatMult_SeqAIJ, 1732 MatMultAdd_SeqAIJ, 1733 MatMultTrans_SeqAIJ, 1734 MatMultTransAdd_SeqAIJ, 1735 MatSolve_SeqAIJ, 1736 MatSolveAdd_SeqAIJ, 1737 MatSolveTrans_SeqAIJ, 1738 MatSolveTransAdd_SeqAIJ, 1739 MatLUFactor_SeqAIJ, 1740 0, 1741 MatRelax_SeqAIJ, 1742 MatTranspose_SeqAIJ, 1743 MatGetInfo_SeqAIJ, 1744 MatEqual_SeqAIJ, 1745 MatGetDiagonal_SeqAIJ, 1746 MatDiagonalScale_SeqAIJ, 1747 MatNorm_SeqAIJ, 1748 0, 1749 MatAssemblyEnd_SeqAIJ, 1750 MatCompress_SeqAIJ, 1751 MatSetOption_SeqAIJ, 1752 MatZeroEntries_SeqAIJ, 1753 MatZeroRows_SeqAIJ, 1754 MatLUFactorSymbolic_SeqAIJ, 1755 MatLUFactorNumeric_SeqAIJ, 1756 0, 1757 0, 1758 MatGetSize_SeqAIJ, 1759 MatGetSize_SeqAIJ, 1760 MatGetOwnershipRange_SeqAIJ, 1761 MatILUFactorSymbolic_SeqAIJ, 1762 0, 1763 0, 1764 0, 1765 MatDuplicate_SeqAIJ, 1766 0, 1767 0, 1768 MatILUFactor_SeqAIJ, 1769 0, 1770 0, 1771 MatGetSubMatrices_SeqAIJ, 1772 MatIncreaseOverlap_SeqAIJ, 1773 MatGetValues_SeqAIJ, 1774 MatCopy_SeqAIJ, 1775 MatPrintHelp_SeqAIJ, 1776 MatScale_SeqAIJ, 1777 0, 1778 0, 1779 MatILUDTFactor_SeqAIJ, 1780 MatGetBlockSize_SeqAIJ, 1781 MatGetRowIJ_SeqAIJ, 1782 MatRestoreRowIJ_SeqAIJ, 1783 MatGetColumnIJ_SeqAIJ, 1784 MatRestoreColumnIJ_SeqAIJ, 1785 MatFDColoringCreate_SeqAIJ, 1786 MatColoringPatch_SeqAIJ, 1787 0, 1788 MatPermute_SeqAIJ, 1789 0, 1790 0, 1791 0, 1792 0, 1793 MatGetMaps_Petsc}; 1794 1795 extern int MatUseSuperLU_SeqAIJ(Mat); 1796 extern int MatUseEssl_SeqAIJ(Mat); 1797 extern int MatUseDXML_SeqAIJ(Mat); 1798 1799 #undef __FUNC__ 1800 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1801 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1802 { 1803 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1804 int i,nz,n; 1805 1806 PetscFunctionBegin; 1807 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1808 1809 nz = aij->maxnz; 1810 n = aij->n; 1811 for (i=0; i<nz; i++) { 1812 aij->j[i] = indices[i]; 1813 } 1814 aij->nz = nz; 1815 for ( i=0; i<n; i++ ) { 1816 aij->ilen[i] = aij->imax[i]; 1817 } 1818 1819 PetscFunctionReturn(0); 1820 } 1821 1822 #undef __FUNC__ 1823 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1824 /*@ 1825 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1826 in the matrix. 1827 1828 Input Parameters: 1829 + mat - the SeqAIJ matrix 1830 - indices - the column indices 1831 1832 Notes: 1833 This can be called if you have precomputed the nonzero structure of the 1834 matrix and want to provide it to the matrix object to improve the performance 1835 of the MatSetValues() operation. 1836 1837 You MUST have set the correct numbers of nonzeros per row in the call to 1838 MatCreateSeqAIJ(). 1839 1840 MUST be called before any calls to MatSetValues(); 1841 1842 @*/ 1843 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1844 { 1845 int ierr,(*f)(Mat,int *); 1846 1847 PetscFunctionBegin; 1848 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1849 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1850 CHKERRQ(ierr); 1851 if (f) { 1852 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1853 } else { 1854 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1855 } 1856 PetscFunctionReturn(0); 1857 } 1858 1859 /* ----------------------------------------------------------------------------------------*/ 1860 1861 #undef __FUNC__ 1862 #define __FUNC__ "MatStoreValues_SeqAIJ" 1863 int MatStoreValues_SeqAIJ(Mat mat) 1864 { 1865 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1866 int nz = aij->i[aij->m]+aij->indexshift; 1867 1868 PetscFunctionBegin; 1869 if (aij->nonew != 1) { 1870 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1871 } 1872 1873 /* allocate space for values if not already there */ 1874 if (!aij->saved_values) { 1875 aij->saved_values = (Scalar *) PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1876 } 1877 1878 /* copy values over */ 1879 PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar)); 1880 PetscFunctionReturn(0); 1881 } 1882 1883 #undef __FUNC__ 1884 #define __FUNC__ "MatStoreValues" 1885 /*@ 1886 MatStoreValues - Stashes a copy of the matrix values; this allows, for 1887 example, reuse of the linear part of a Jacobian, while recomputing the 1888 nonlinear portion. 1889 1890 Collect on Mat 1891 1892 Input Parameters: 1893 . mat - the matrix (currently on AIJ matrices support this option) 1894 1895 Common Usage, with SNESSolve(): 1896 $ Create Jacobian matrix 1897 $ Set linear terms into matrix 1898 $ Apply boundary conditions to matrix, at this time matrix must have 1899 $ final nonzero structure (i.e. setting the nonlinear terms and applying 1900 $ boundary conditions again will not change the nonzero structure 1901 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1902 $ ierr = MatStoreValues(mat); 1903 $ Call SNESSetJacobian() with matrix 1904 $ In your Jacobian routine 1905 $ ierr = MatRetrieveValues(mat); 1906 $ Set nonlinear terms in matrix 1907 1908 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 1909 $ // build linear portion of Jacobian 1910 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1911 $ ierr = MatStoreValues(mat); 1912 $ loop over nonlinear iterations 1913 $ ierr = MatRetrieveValues(mat); 1914 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 1915 $ // call MatAssemblyBegin/End() on matrix 1916 $ Solve linear system with Jacobian 1917 $ endloop 1918 1919 Notes: 1920 Matrix must already be assemblied before calling this routine 1921 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 1922 calling this routine. 1923 1924 .seealso: MatRetrieveValues() 1925 1926 @*/ 1927 int MatStoreValues(Mat mat) 1928 { 1929 int ierr,(*f)(Mat); 1930 1931 PetscFunctionBegin; 1932 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1933 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1934 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1935 1936 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f); 1937 CHKERRQ(ierr); 1938 if (f) { 1939 ierr = (*f)(mat);CHKERRQ(ierr); 1940 } else { 1941 SETERRQ(1,1,"Wrong type of matrix to store values"); 1942 } 1943 PetscFunctionReturn(0); 1944 } 1945 1946 #undef __FUNC__ 1947 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 1948 int MatRetrieveValues_SeqAIJ(Mat mat) 1949 { 1950 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1951 int nz = aij->i[aij->m]+aij->indexshift; 1952 1953 PetscFunctionBegin; 1954 if (aij->nonew != 1) { 1955 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1956 } 1957 if (!aij->saved_values) { 1958 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1959 } 1960 1961 /* copy values over */ 1962 PetscMemcpy(aij->a, aij->saved_values,nz*sizeof(Scalar)); 1963 PetscFunctionReturn(0); 1964 } 1965 1966 #undef __FUNC__ 1967 #define __FUNC__ "MatRetrieveValues" 1968 /*@ 1969 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 1970 example, reuse of the linear part of a Jacobian, while recomputing the 1971 nonlinear portion. 1972 1973 Collect on Mat 1974 1975 Input Parameters: 1976 . mat - the matrix (currently on AIJ matrices support this option) 1977 1978 .seealso: MatStoreValues() 1979 1980 @*/ 1981 int MatRetrieveValues(Mat mat) 1982 { 1983 int ierr,(*f)(Mat); 1984 1985 PetscFunctionBegin; 1986 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1987 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1988 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1989 1990 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f); 1991 CHKERRQ(ierr); 1992 if (f) { 1993 ierr = (*f)(mat);CHKERRQ(ierr); 1994 } else { 1995 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 1996 } 1997 PetscFunctionReturn(0); 1998 } 1999 2000 /* --------------------------------------------------------------------------------*/ 2001 2002 #undef __FUNC__ 2003 #define __FUNC__ "MatCreateSeqAIJ" 2004 /*@C 2005 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2006 (the default parallel PETSc format). For good matrix assembly performance 2007 the user should preallocate the matrix storage by setting the parameter nz 2008 (or the array nzz). By setting these parameters accurately, performance 2009 during matrix assembly can be increased by more than a factor of 50. 2010 2011 Collective on MPI_Comm 2012 2013 Input Parameters: 2014 + comm - MPI communicator, set to PETSC_COMM_SELF 2015 . m - number of rows 2016 . n - number of columns 2017 . nz - number of nonzeros per row (same for all rows) 2018 - nzz - array containing the number of nonzeros in the various rows 2019 (possibly different for each row) or PETSC_NULL 2020 2021 Output Parameter: 2022 . A - the matrix 2023 2024 Notes: 2025 The AIJ format (also called the Yale sparse matrix format or 2026 compressed row storage), is fully compatible with standard Fortran 77 2027 storage. That is, the stored row and column indices can begin at 2028 either one (as in Fortran) or zero. See the users' manual for details. 2029 2030 Specify the preallocated storage with either nz or nnz (not both). 2031 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2032 allocation. For large problems you MUST preallocate memory or you 2033 will get TERRIBLE performance, see the users' manual chapter on matrices. 2034 2035 By default, this format uses inodes (identical nodes) when possible, to 2036 improve numerical efficiency of matrix-vector products and solves. We 2037 search for consecutive rows with the same nonzero structure, thereby 2038 reusing matrix information to achieve increased efficiency. 2039 2040 Options Database Keys: 2041 + -mat_aij_no_inode - Do not use inodes 2042 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2043 - -mat_aij_oneindex - Internally use indexing starting at 1 2044 rather than 0. Note that when calling MatSetValues(), 2045 the user still MUST index entries starting at 0! 2046 2047 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 2048 @*/ 2049 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 2050 { 2051 Mat B; 2052 Mat_SeqAIJ *b; 2053 int i, len, ierr, flg,size; 2054 2055 PetscFunctionBegin; 2056 MPI_Comm_size(comm,&size); 2057 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2058 2059 *A = 0; 2060 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 2061 PLogObjectCreate(B); 2062 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 2063 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 2064 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 2065 B->ops->destroy = MatDestroy_SeqAIJ; 2066 B->ops->view = MatView_SeqAIJ; 2067 B->factor = 0; 2068 B->lupivotthreshold = 1.0; 2069 B->mapping = 0; 2070 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 2071 b->ilu_preserve_row_sums = PETSC_FALSE; 2072 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2073 b->row = 0; 2074 b->col = 0; 2075 b->icol = 0; 2076 b->indexshift = 0; 2077 b->reallocs = 0; 2078 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 2079 if (flg) b->indexshift = -1; 2080 2081 b->m = m; B->m = m; B->M = m; 2082 b->n = n; B->n = n; B->N = n; 2083 2084 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2085 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2086 2087 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 2088 if (nnz == PETSC_NULL) { 2089 if (nz == PETSC_DEFAULT) nz = 10; 2090 else if (nz <= 0) nz = 1; 2091 for ( i=0; i<m; i++ ) b->imax[i] = nz; 2092 nz = nz*m; 2093 } else { 2094 nz = 0; 2095 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 2096 } 2097 2098 /* allocate the matrix space */ 2099 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2100 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 2101 b->j = (int *) (b->a + nz); 2102 PetscMemzero(b->j,nz*sizeof(int)); 2103 b->i = b->j + nz; 2104 b->singlemalloc = 1; 2105 2106 b->i[0] = -b->indexshift; 2107 for (i=1; i<m+1; i++) { 2108 b->i[i] = b->i[i-1] + b->imax[i-1]; 2109 } 2110 2111 /* b->ilen will count nonzeros in each row so far. */ 2112 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 2113 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2114 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 2115 2116 b->nz = 0; 2117 b->maxnz = nz; 2118 b->sorted = 0; 2119 b->roworiented = 1; 2120 b->nonew = 0; 2121 b->diag = 0; 2122 b->solve_work = 0; 2123 b->spptr = 0; 2124 b->inode.node_count = 0; 2125 b->inode.size = 0; 2126 b->inode.limit = 5; 2127 b->inode.max_limit = 5; 2128 b->saved_values = 0; 2129 B->info.nz_unneeded = (double)b->maxnz; 2130 2131 *A = B; 2132 2133 /* SuperLU is not currently supported through PETSc */ 2134 #if defined(HAVE_SUPERLU) 2135 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 2136 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 2137 #endif 2138 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 2139 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 2140 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 2141 if (flg) { 2142 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2143 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 2144 } 2145 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 2146 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 2147 2148 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2149 "MatSeqAIJSetColumnIndices_SeqAIJ", 2150 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2151 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C", 2152 "MatStoreValues_SeqAIJ", 2153 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2154 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C", 2155 "MatRetrieveValues_SeqAIJ", 2156 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2157 PetscFunctionReturn(0); 2158 } 2159 2160 #undef __FUNC__ 2161 #define __FUNC__ "MatDuplicate_SeqAIJ" 2162 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2163 { 2164 Mat C; 2165 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 2166 int i,len, m = a->m,shift = a->indexshift,ierr; 2167 2168 PetscFunctionBegin; 2169 *B = 0; 2170 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 2171 PLogObjectCreate(C); 2172 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 2173 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 2174 C->ops->destroy = MatDestroy_SeqAIJ; 2175 C->ops->view = MatView_SeqAIJ; 2176 C->factor = A->factor; 2177 c->row = 0; 2178 c->col = 0; 2179 c->icol = 0; 2180 c->indexshift = shift; 2181 C->assembled = PETSC_TRUE; 2182 2183 c->m = C->m = a->m; 2184 c->n = C->n = a->n; 2185 C->M = a->m; 2186 C->N = a->n; 2187 2188 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 2189 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 2190 for ( i=0; i<m; i++ ) { 2191 c->imax[i] = a->imax[i]; 2192 c->ilen[i] = a->ilen[i]; 2193 } 2194 2195 /* allocate the matrix space */ 2196 c->singlemalloc = 1; 2197 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2198 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 2199 c->j = (int *) (c->a + a->i[m] + shift); 2200 c->i = c->j + a->i[m] + shift; 2201 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 2202 if (m > 0) { 2203 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 2204 if (cpvalues == MAT_COPY_VALUES) { 2205 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 2206 } else { 2207 PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar)); 2208 } 2209 } 2210 2211 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2212 c->sorted = a->sorted; 2213 c->roworiented = a->roworiented; 2214 c->nonew = a->nonew; 2215 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2216 c->saved_values = 0; 2217 2218 if (a->diag) { 2219 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2220 PLogObjectMemory(C,(m+1)*sizeof(int)); 2221 for ( i=0; i<m; i++ ) { 2222 c->diag[i] = a->diag[i]; 2223 } 2224 } else c->diag = 0; 2225 c->inode.limit = a->inode.limit; 2226 c->inode.max_limit = a->inode.max_limit; 2227 if (a->inode.size){ 2228 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2229 c->inode.node_count = a->inode.node_count; 2230 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2231 } else { 2232 c->inode.size = 0; 2233 c->inode.node_count = 0; 2234 } 2235 c->nz = a->nz; 2236 c->maxnz = a->maxnz; 2237 c->solve_work = 0; 2238 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2239 2240 *B = C; 2241 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2242 "MatSeqAIJSetColumnIndices_SeqAIJ", 2243 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2244 PetscFunctionReturn(0); 2245 } 2246 2247 #undef __FUNC__ 2248 #define __FUNC__ "MatLoad_SeqAIJ" 2249 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2250 { 2251 Mat_SeqAIJ *a; 2252 Mat B; 2253 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2254 MPI_Comm comm; 2255 2256 PetscFunctionBegin; 2257 PetscObjectGetComm((PetscObject) viewer,&comm); 2258 MPI_Comm_size(comm,&size); 2259 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2260 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2261 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2262 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2263 M = header[1]; N = header[2]; nz = header[3]; 2264 2265 if (nz < 0) { 2266 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2267 } 2268 2269 /* read in row lengths */ 2270 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2271 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2272 2273 /* create our matrix */ 2274 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2275 B = *A; 2276 a = (Mat_SeqAIJ *) B->data; 2277 shift = a->indexshift; 2278 2279 /* read in column indices and adjust for Fortran indexing*/ 2280 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2281 if (shift) { 2282 for ( i=0; i<nz; i++ ) { 2283 a->j[i] += 1; 2284 } 2285 } 2286 2287 /* read in nonzero values */ 2288 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2289 2290 /* set matrix "i" values */ 2291 a->i[0] = -shift; 2292 for ( i=1; i<= M; i++ ) { 2293 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2294 a->ilen[i-1] = rowlengths[i-1]; 2295 } 2296 PetscFree(rowlengths); 2297 2298 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2299 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2300 PetscFunctionReturn(0); 2301 } 2302 2303 #undef __FUNC__ 2304 #define __FUNC__ "MatEqual_SeqAIJ" 2305 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2306 { 2307 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2308 2309 PetscFunctionBegin; 2310 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2311 2312 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2313 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2314 (a->indexshift != b->indexshift)) { 2315 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2316 } 2317 2318 /* if the a->i are the same */ 2319 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2320 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2321 } 2322 2323 /* if a->j are the same */ 2324 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2325 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2326 } 2327 2328 /* if a->a are the same */ 2329 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2330 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2331 } 2332 *flg = PETSC_TRUE; 2333 PetscFunctionReturn(0); 2334 2335 } 2336