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