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