1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.278 1998/07/14 21:27:22 bsmith Exp bsmith $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "src/inline/dot.h" 16 #include "src/inline/bitarray.h" 17 18 /* 19 Basic AIJ format ILU based on drop tolerance 20 */ 21 #undef __FUNC__ 22 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 23 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 24 { 25 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 26 int ierr = 1; 27 28 PetscFunctionBegin; 29 SETERRQ(ierr,0,"Not implemented"); 30 #if !defined(USE_PETSC_DEBUG) 31 PetscFunctionReturn(0); 32 #endif 33 } 34 35 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 36 37 #undef __FUNC__ 38 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 39 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 40 PetscTruth *done) 41 { 42 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 43 int ierr,i,ishift; 44 45 PetscFunctionBegin; 46 *m = A->m; 47 if (!ia) PetscFunctionReturn(0); 48 ishift = a->indexshift; 49 if (symmetric) { 50 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 51 } else if (oshift == 0 && ishift == -1) { 52 int nz = a->i[a->m]; 53 /* malloc space and subtract 1 from i and j indices */ 54 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 55 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 56 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 57 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 58 } else if (oshift == 1 && ishift == 0) { 59 int nz = a->i[a->m] + 1; 60 /* malloc space and add 1 to i and j indices */ 61 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 62 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 63 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 64 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 65 } else { 66 *ia = a->i; *ja = a->j; 67 } 68 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNC__ 73 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 74 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 75 PetscTruth *done) 76 { 77 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 78 int ishift = a->indexshift; 79 80 PetscFunctionBegin; 81 if (!ia) PetscFunctionReturn(0); 82 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 83 PetscFree(*ia); 84 PetscFree(*ja); 85 } 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNC__ 90 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 91 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 92 PetscTruth *done) 93 { 94 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 95 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 96 int nz = a->i[m]+ishift,row,*jj,mr,col; 97 98 PetscFunctionBegin; 99 *nn = A->n; 100 if (!ia) PetscFunctionReturn(0); 101 if (symmetric) { 102 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 103 } else { 104 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 105 PetscMemzero(collengths,n*sizeof(int)); 106 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 107 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 108 jj = a->j; 109 for ( i=0; i<nz; i++ ) { 110 collengths[jj[i] + ishift]++; 111 } 112 cia[0] = oshift; 113 for ( i=0; i<n; i++) { 114 cia[i+1] = cia[i] + collengths[i]; 115 } 116 PetscMemzero(collengths,n*sizeof(int)); 117 jj = a->j; 118 for ( row=0; row<m; row++ ) { 119 mr = a->i[row+1] - a->i[row]; 120 for ( i=0; i<mr; i++ ) { 121 col = *jj++ + ishift; 122 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 123 } 124 } 125 PetscFree(collengths); 126 *ia = cia; *ja = cja; 127 } 128 129 PetscFunctionReturn(0); 130 } 131 132 #undef __FUNC__ 133 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 134 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 135 int **ja,PetscTruth *done) 136 { 137 PetscFunctionBegin; 138 if (!ia) PetscFunctionReturn(0); 139 140 PetscFree(*ia); 141 PetscFree(*ja); 142 143 PetscFunctionReturn(0); 144 } 145 146 #define CHUNKSIZE 15 147 148 #undef __FUNC__ 149 #define __FUNC__ "MatSetValues_SeqAIJ" 150 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 151 { 152 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 153 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 154 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 155 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 156 Scalar *ap,value, *aa = a->a; 157 158 PetscFunctionBegin; 159 for ( k=0; k<m; k++ ) { /* loop over added rows */ 160 row = im[k]; 161 #if defined(USE_PETSC_BOPT_g) 162 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 163 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 164 #endif 165 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 166 rmax = imax[row]; nrow = ailen[row]; 167 low = 0; 168 for ( l=0; l<n; l++ ) { /* loop over added columns */ 169 #if defined(USE_PETSC_BOPT_g) 170 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 171 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 172 #endif 173 col = in[l] - shift; 174 if (roworiented) { 175 value = *v++; 176 } else { 177 value = v[k + l*m]; 178 } 179 if (!sorted) low = 0; high = nrow; 180 while (high-low > 5) { 181 t = (low+high)/2; 182 if (rp[t] > col) high = t; 183 else low = t; 184 } 185 for ( i=low; i<high; i++ ) { 186 if (rp[i] > col) break; 187 if (rp[i] == col) { 188 if (is == ADD_VALUES) ap[i] += value; 189 else ap[i] = value; 190 goto noinsert; 191 } 192 } 193 if (nonew == 1) goto noinsert; 194 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 195 if (nrow >= rmax) { 196 /* there is no extra room in row, therefore enlarge */ 197 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 198 Scalar *new_a; 199 200 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 201 202 /* malloc new storage space */ 203 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 204 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 205 new_j = (int *) (new_a + new_nz); 206 new_i = new_j + new_nz; 207 208 /* copy over old data into new slots */ 209 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 210 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 211 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 212 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 213 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 214 len*sizeof(int)); 215 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 216 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 217 len*sizeof(Scalar)); 218 /* free up old matrix storage */ 219 PetscFree(a->a); 220 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 221 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 222 a->singlemalloc = 1; 223 224 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 225 rmax = imax[row] = imax[row] + CHUNKSIZE; 226 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 227 a->maxnz += CHUNKSIZE; 228 a->reallocs++; 229 } 230 N = nrow++ - 1; a->nz++; 231 /* shift up all the later entries in this row */ 232 for ( ii=N; ii>=i; ii-- ) { 233 rp[ii+1] = rp[ii]; 234 ap[ii+1] = ap[ii]; 235 } 236 rp[i] = col; 237 ap[i] = value; 238 noinsert:; 239 low = i + 1; 240 } 241 ailen[row] = nrow; 242 } 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNC__ 247 #define __FUNC__ "MatGetValues_SeqAIJ" 248 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 249 { 250 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 251 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 252 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 253 Scalar *ap, *aa = a->a, zero = 0.0; 254 255 PetscFunctionBegin; 256 for ( k=0; k<m; k++ ) { /* loop over rows */ 257 row = im[k]; 258 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 259 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 260 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 261 nrow = ailen[row]; 262 for ( l=0; l<n; l++ ) { /* loop over columns */ 263 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 264 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 265 col = in[l] - shift; 266 high = nrow; low = 0; /* assume unsorted */ 267 while (high-low > 5) { 268 t = (low+high)/2; 269 if (rp[t] > col) high = t; 270 else low = t; 271 } 272 for ( i=low; i<high; i++ ) { 273 if (rp[i] > col) break; 274 if (rp[i] == col) { 275 *v++ = ap[i]; 276 goto finished; 277 } 278 } 279 *v++ = zero; 280 finished:; 281 } 282 } 283 PetscFunctionReturn(0); 284 } 285 286 287 #undef __FUNC__ 288 #define __FUNC__ "MatView_SeqAIJ_Binary" 289 int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 290 { 291 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 292 int i, fd, *col_lens, ierr; 293 294 PetscFunctionBegin; 295 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 296 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 297 col_lens[0] = MAT_COOKIE; 298 col_lens[1] = a->m; 299 col_lens[2] = a->n; 300 col_lens[3] = a->nz; 301 302 /* store lengths of each row and write (including header) to file */ 303 for ( i=0; i<a->m; i++ ) { 304 col_lens[4+i] = a->i[i+1] - a->i[i]; 305 } 306 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 307 PetscFree(col_lens); 308 309 /* store column indices (zero start index) */ 310 if (a->indexshift) { 311 for ( i=0; i<a->nz; i++ ) a->j[i]--; 312 } 313 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0); CHKERRQ(ierr); 314 if (a->indexshift) { 315 for ( i=0; i<a->nz; i++ ) a->j[i]++; 316 } 317 318 /* store nonzero values */ 319 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 #undef __FUNC__ 324 #define __FUNC__ "MatView_SeqAIJ_ASCII" 325 int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 326 { 327 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 328 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 329 FILE *fd; 330 char *outputname; 331 332 PetscFunctionBegin; 333 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 334 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 335 ierr = ViewerGetFormat(viewer,&format); 336 if (format == VIEWER_FORMAT_ASCII_INFO) { 337 PetscFunctionReturn(0); 338 } else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 339 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 340 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 341 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 342 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 343 a->inode.node_count,a->inode.limit); 344 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 345 int nofinalvalue = 0; 346 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != a->n-!shift)) { 347 nofinalvalue = 1; 348 } 349 fprintf(fd,"%% Size = %d %d \n",m,a->n); 350 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 351 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue); 352 fprintf(fd,"zzz = [\n"); 353 354 for (i=0; i<m; i++) { 355 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 356 #if defined(USE_PETSC_COMPLEX) 357 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,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 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 916 if (a->j[j]+shift == i) { 917 diag[i] = j - shift; 918 break; 919 } 920 } 921 } 922 a->diag = diag; 923 PetscFunctionReturn(0); 924 } 925 926 #undef __FUNC__ 927 #define __FUNC__ "MatRelax_SeqAIJ" 928 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 929 double fshift,int its,Vec xx) 930 { 931 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 932 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 933 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 934 935 PetscFunctionBegin; 936 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 937 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 938 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 939 diag = a->diag; 940 xs = x + shift; /* shifted by one for index start of a or a->j*/ 941 if (flag == SOR_APPLY_UPPER) { 942 /* apply ( U + D/omega) to the vector */ 943 bs = b + shift; 944 for ( i=0; i<m; i++ ) { 945 d = fshift + a->a[diag[i] + shift]; 946 n = a->i[i+1] - diag[i] - 1; 947 PLogFlops(2*n-1); 948 idx = a->j + diag[i] + (!shift); 949 v = a->a + diag[i] + (!shift); 950 sum = b[i]*d/omega; 951 SPARSEDENSEDOT(sum,bs,v,idx,n); 952 x[i] = sum; 953 } 954 PetscFunctionReturn(0); 955 } 956 if (flag == SOR_APPLY_LOWER) { 957 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 958 } else if (flag & SOR_EISENSTAT) { 959 /* Let A = L + U + D; where L is lower trianglar, 960 U is upper triangular, E is diagonal; This routine applies 961 962 (L + E)^{-1} A (U + E)^{-1} 963 964 to a vector efficiently using Eisenstat's trick. This is for 965 the case of SSOR preconditioner, so E is D/omega where omega 966 is the relaxation factor. 967 */ 968 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 969 scale = (2.0/omega) - 1.0; 970 971 /* x = (E + U)^{-1} b */ 972 for ( i=m-1; i>=0; i-- ) { 973 d = fshift + a->a[diag[i] + shift]; 974 n = a->i[i+1] - diag[i] - 1; 975 PLogFlops(2*n-1); 976 idx = a->j + diag[i] + (!shift); 977 v = a->a + diag[i] + (!shift); 978 sum = b[i]; 979 SPARSEDENSEMDOT(sum,xs,v,idx,n); 980 x[i] = omega*(sum/d); 981 } 982 983 /* t = b - (2*E - D)x */ 984 v = a->a; 985 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 986 PLogFlops(2*m); 987 988 989 /* t = (E + L)^{-1}t */ 990 ts = t + shift; /* shifted by one for index start of a or a->j*/ 991 diag = a->diag; 992 for ( i=0; i<m; i++ ) { 993 d = fshift + a->a[diag[i]+shift]; 994 n = diag[i] - a->i[i]; 995 PLogFlops(2*n-1); 996 idx = a->j + a->i[i] + shift; 997 v = a->a + a->i[i] + shift; 998 sum = t[i]; 999 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1000 t[i] = omega*(sum/d); 1001 } 1002 1003 /* x = x + t */ 1004 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1005 PetscFree(t); 1006 PetscFunctionReturn(0); 1007 } 1008 if (flag & SOR_ZERO_INITIAL_GUESS) { 1009 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1010 for ( i=0; i<m; i++ ) { 1011 d = fshift + a->a[diag[i]+shift]; 1012 n = diag[i] - a->i[i]; 1013 PLogFlops(2*n-1); 1014 idx = a->j + a->i[i] + shift; 1015 v = a->a + a->i[i] + shift; 1016 sum = b[i]; 1017 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1018 x[i] = omega*(sum/d); 1019 } 1020 xb = x; 1021 } else xb = b; 1022 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1023 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1024 for ( i=0; i<m; i++ ) { 1025 x[i] *= a->a[diag[i]+shift]; 1026 } 1027 PLogFlops(m); 1028 } 1029 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1030 for ( i=m-1; i>=0; i-- ) { 1031 d = fshift + a->a[diag[i] + shift]; 1032 n = a->i[i+1] - diag[i] - 1; 1033 PLogFlops(2*n-1); 1034 idx = a->j + diag[i] + (!shift); 1035 v = a->a + diag[i] + (!shift); 1036 sum = xb[i]; 1037 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1038 x[i] = omega*(sum/d); 1039 } 1040 } 1041 its--; 1042 } 1043 while (its--) { 1044 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1045 for ( i=0; i<m; i++ ) { 1046 d = fshift + a->a[diag[i]+shift]; 1047 n = a->i[i+1] - a->i[i]; 1048 PLogFlops(2*n-1); 1049 idx = a->j + a->i[i] + shift; 1050 v = a->a + a->i[i] + shift; 1051 sum = b[i]; 1052 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1053 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1054 } 1055 } 1056 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1057 for ( i=m-1; i>=0; i-- ) { 1058 d = fshift + a->a[diag[i] + shift]; 1059 n = a->i[i+1] - a->i[i]; 1060 PLogFlops(2*n-1); 1061 idx = a->j + a->i[i] + shift; 1062 v = a->a + a->i[i] + shift; 1063 sum = b[i]; 1064 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1065 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1066 } 1067 } 1068 } 1069 PetscFunctionReturn(0); 1070 } 1071 1072 #undef __FUNC__ 1073 #define __FUNC__ "MatGetInfo_SeqAIJ" 1074 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1075 { 1076 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1077 1078 PetscFunctionBegin; 1079 info->rows_global = (double)a->m; 1080 info->columns_global = (double)a->n; 1081 info->rows_local = (double)a->m; 1082 info->columns_local = (double)a->n; 1083 info->block_size = 1.0; 1084 info->nz_allocated = (double)a->maxnz; 1085 info->nz_used = (double)a->nz; 1086 info->nz_unneeded = (double)(a->maxnz - a->nz); 1087 info->assemblies = (double)A->num_ass; 1088 info->mallocs = (double)a->reallocs; 1089 info->memory = A->mem; 1090 if (A->factor) { 1091 info->fill_ratio_given = A->info.fill_ratio_given; 1092 info->fill_ratio_needed = A->info.fill_ratio_needed; 1093 info->factor_mallocs = A->info.factor_mallocs; 1094 } else { 1095 info->fill_ratio_given = 0; 1096 info->fill_ratio_needed = 0; 1097 info->factor_mallocs = 0; 1098 } 1099 PetscFunctionReturn(0); 1100 } 1101 1102 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1103 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1104 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1105 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1106 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1107 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1108 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1109 1110 #undef __FUNC__ 1111 #define __FUNC__ "MatZeroRows_SeqAIJ" 1112 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1113 { 1114 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1115 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1116 1117 PetscFunctionBegin; 1118 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1119 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1120 if (diag) { 1121 for ( i=0; i<N; i++ ) { 1122 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1123 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1124 a->ilen[rows[i]] = 1; 1125 a->a[a->i[rows[i]]+shift] = *diag; 1126 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1127 } else { 1128 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1129 } 1130 } 1131 } else { 1132 for ( i=0; i<N; i++ ) { 1133 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1134 a->ilen[rows[i]] = 0; 1135 } 1136 } 1137 ISRestoreIndices(is,&rows); 1138 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #undef __FUNC__ 1143 #define __FUNC__ "MatGetSize_SeqAIJ" 1144 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1145 { 1146 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1147 1148 PetscFunctionBegin; 1149 if (m) *m = a->m; 1150 if (n) *n = a->n; 1151 PetscFunctionReturn(0); 1152 } 1153 1154 #undef __FUNC__ 1155 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1156 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1157 { 1158 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1159 1160 PetscFunctionBegin; 1161 *m = 0; *n = a->m; 1162 PetscFunctionReturn(0); 1163 } 1164 1165 #undef __FUNC__ 1166 #define __FUNC__ "MatGetRow_SeqAIJ" 1167 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1168 { 1169 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1170 int *itmp,i,shift = a->indexshift; 1171 1172 PetscFunctionBegin; 1173 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1174 1175 *nz = a->i[row+1] - a->i[row]; 1176 if (v) *v = a->a + a->i[row] + shift; 1177 if (idx) { 1178 itmp = a->j + a->i[row] + shift; 1179 if (*nz && shift) { 1180 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1181 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1182 } else if (*nz) { 1183 *idx = itmp; 1184 } 1185 else *idx = 0; 1186 } 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNC__ 1191 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1192 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1193 { 1194 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1195 1196 PetscFunctionBegin; 1197 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1198 PetscFunctionReturn(0); 1199 } 1200 1201 #undef __FUNC__ 1202 #define __FUNC__ "MatNorm_SeqAIJ" 1203 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1204 { 1205 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1206 Scalar *v = a->a; 1207 double sum = 0.0; 1208 int i, j,shift = a->indexshift; 1209 1210 PetscFunctionBegin; 1211 if (type == NORM_FROBENIUS) { 1212 for (i=0; i<a->nz; i++ ) { 1213 #if defined(USE_PETSC_COMPLEX) 1214 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1215 #else 1216 sum += (*v)*(*v); v++; 1217 #endif 1218 } 1219 *norm = sqrt(sum); 1220 } else if (type == NORM_1) { 1221 double *tmp; 1222 int *jj = a->j; 1223 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1224 PetscMemzero(tmp,a->n*sizeof(double)); 1225 *norm = 0.0; 1226 for ( j=0; j<a->nz; j++ ) { 1227 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1228 } 1229 for ( j=0; j<a->n; j++ ) { 1230 if (tmp[j] > *norm) *norm = tmp[j]; 1231 } 1232 PetscFree(tmp); 1233 } else if (type == NORM_INFINITY) { 1234 *norm = 0.0; 1235 for ( j=0; j<a->m; j++ ) { 1236 v = a->a + a->i[j] + shift; 1237 sum = 0.0; 1238 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1239 sum += PetscAbsScalar(*v); v++; 1240 } 1241 if (sum > *norm) *norm = sum; 1242 } 1243 } else { 1244 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1245 } 1246 PetscFunctionReturn(0); 1247 } 1248 1249 #undef __FUNC__ 1250 #define __FUNC__ "MatTranspose_SeqAIJ" 1251 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1252 { 1253 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1254 Mat C; 1255 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1256 int shift = a->indexshift; 1257 Scalar *array = a->a; 1258 1259 PetscFunctionBegin; 1260 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1261 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1262 PetscMemzero(col,(1+a->n)*sizeof(int)); 1263 if (shift) { 1264 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1265 } 1266 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1267 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1268 PetscFree(col); 1269 for ( i=0; i<m; i++ ) { 1270 len = ai[i+1]-ai[i]; 1271 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1272 array += len; aj += len; 1273 } 1274 if (shift) { 1275 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1276 } 1277 1278 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1279 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1280 1281 if (B != PETSC_NULL) { 1282 *B = C; 1283 } else { 1284 PetscOps *Abops; 1285 MatOps Aops; 1286 1287 /* This isn't really an in-place transpose */ 1288 PetscFree(a->a); 1289 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1290 if (a->diag) PetscFree(a->diag); 1291 if (a->ilen) PetscFree(a->ilen); 1292 if (a->imax) PetscFree(a->imax); 1293 if (a->solve_work) PetscFree(a->solve_work); 1294 if (a->inode.size) PetscFree(a->inode.size); 1295 PetscFree(a); 1296 1297 1298 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1299 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1300 1301 /* 1302 This is horrible, horrible code. We need to keep the 1303 the bops and ops Structures, copy everything from C 1304 including the function pointers.. 1305 */ 1306 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1307 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1308 Abops = A->bops; 1309 Aops = A->ops; 1310 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1311 A->bops = Abops; 1312 A->ops = Aops; 1313 A->qlist = 0; 1314 1315 PetscHeaderDestroy(C); 1316 } 1317 PetscFunctionReturn(0); 1318 } 1319 1320 #undef __FUNC__ 1321 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1322 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1323 { 1324 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1325 Scalar *l,*r,x,*v; 1326 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1327 1328 PetscFunctionBegin; 1329 if (ll) { 1330 /* The local size is used so that VecMPI can be passed to this routine 1331 by MatDiagonalScale_MPIAIJ */ 1332 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1333 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1334 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1335 v = a->a; 1336 for ( i=0; i<m; i++ ) { 1337 x = l[i]; 1338 M = a->i[i+1] - a->i[i]; 1339 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1340 } 1341 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1342 PLogFlops(nz); 1343 } 1344 if (rr) { 1345 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1346 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1347 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1348 v = a->a; jj = a->j; 1349 for ( i=0; i<nz; i++ ) { 1350 (*v++) *= r[*jj++ + shift]; 1351 } 1352 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1353 PLogFlops(nz); 1354 } 1355 PetscFunctionReturn(0); 1356 } 1357 1358 #undef __FUNC__ 1359 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1360 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1361 { 1362 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1363 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1364 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1365 register int sum,lensi; 1366 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1367 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1368 Scalar *a_new,*mat_a; 1369 Mat C; 1370 PetscTruth stride; 1371 1372 PetscFunctionBegin; 1373 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1374 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1375 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1376 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1377 1378 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1379 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1380 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1381 1382 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1383 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1384 if (stride && step == 1) { 1385 /* special case of contiguous rows */ 1386 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1387 starts = lens + ncols; 1388 /* loop over new rows determining lens and starting points */ 1389 for (i=0; i<nrows; i++) { 1390 kstart = ai[irow[i]]+shift; 1391 kend = kstart + ailen[irow[i]]; 1392 for ( k=kstart; k<kend; k++ ) { 1393 if (aj[k]+shift >= first) { 1394 starts[i] = k; 1395 break; 1396 } 1397 } 1398 sum = 0; 1399 while (k < kend) { 1400 if (aj[k++]+shift >= first+ncols) break; 1401 sum++; 1402 } 1403 lens[i] = sum; 1404 } 1405 /* create submatrix */ 1406 if (scall == MAT_REUSE_MATRIX) { 1407 int n_cols,n_rows; 1408 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1409 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1410 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1411 C = *B; 1412 } else { 1413 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1414 } 1415 c = (Mat_SeqAIJ*) C->data; 1416 1417 /* loop over rows inserting into submatrix */ 1418 a_new = c->a; 1419 j_new = c->j; 1420 i_new = c->i; 1421 i_new[0] = -shift; 1422 for (i=0; i<nrows; i++) { 1423 ii = starts[i]; 1424 lensi = lens[i]; 1425 for ( k=0; k<lensi; k++ ) { 1426 *j_new++ = aj[ii+k] - first; 1427 } 1428 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1429 a_new += lensi; 1430 i_new[i+1] = i_new[i] + lensi; 1431 c->ilen[i] = lensi; 1432 } 1433 PetscFree(lens); 1434 } else { 1435 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1436 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1437 ssmap = smap + shift; 1438 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1439 PetscMemzero(smap,oldcols*sizeof(int)); 1440 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1441 /* determine lens of each row */ 1442 for (i=0; i<nrows; i++) { 1443 kstart = ai[irow[i]]+shift; 1444 kend = kstart + a->ilen[irow[i]]; 1445 lens[i] = 0; 1446 for ( k=kstart; k<kend; k++ ) { 1447 if (ssmap[aj[k]]) { 1448 lens[i]++; 1449 } 1450 } 1451 } 1452 /* Create and fill new matrix */ 1453 if (scall == MAT_REUSE_MATRIX) { 1454 c = (Mat_SeqAIJ *)((*B)->data); 1455 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1456 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1457 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1458 } 1459 PetscMemzero(c->ilen,c->m*sizeof(int)); 1460 C = *B; 1461 } else { 1462 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1463 } 1464 c = (Mat_SeqAIJ *)(C->data); 1465 for (i=0; i<nrows; i++) { 1466 row = irow[i]; 1467 kstart = ai[row]+shift; 1468 kend = kstart + a->ilen[row]; 1469 mat_i = c->i[i]+shift; 1470 mat_j = c->j + mat_i; 1471 mat_a = c->a + mat_i; 1472 mat_ilen = c->ilen + i; 1473 for ( k=kstart; k<kend; k++ ) { 1474 if ((tcol=ssmap[a->j[k]])) { 1475 *mat_j++ = tcol - (!shift); 1476 *mat_a++ = a->a[k]; 1477 (*mat_ilen)++; 1478 1479 } 1480 } 1481 } 1482 /* Free work space */ 1483 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1484 PetscFree(smap); PetscFree(lens); 1485 } 1486 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1487 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1488 1489 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1490 *B = C; 1491 PetscFunctionReturn(0); 1492 } 1493 1494 /* 1495 note: This can only work for identity for row and col. It would 1496 be good to check this and otherwise generate an error. 1497 */ 1498 #undef __FUNC__ 1499 #define __FUNC__ "MatILUFactor_SeqAIJ" 1500 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1501 { 1502 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1503 int ierr; 1504 Mat outA; 1505 1506 PetscFunctionBegin; 1507 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1508 1509 outA = inA; 1510 inA->factor = FACTOR_LU; 1511 a->row = row; 1512 a->col = col; 1513 1514 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1515 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1516 PLogObjectParent(inA,a->icol); 1517 1518 if (!a->solve_work) { /* this matrix may have been factored before */ 1519 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1520 } 1521 1522 if (!a->diag) { 1523 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1524 } 1525 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #include "pinclude/blaslapack.h" 1530 #undef __FUNC__ 1531 #define __FUNC__ "MatScale_SeqAIJ" 1532 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1533 { 1534 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1535 int one = 1; 1536 1537 PetscFunctionBegin; 1538 BLscal_( &a->nz, alpha, a->a, &one ); 1539 PLogFlops(a->nz); 1540 PetscFunctionReturn(0); 1541 } 1542 1543 #undef __FUNC__ 1544 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1545 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1546 { 1547 int ierr,i; 1548 1549 PetscFunctionBegin; 1550 if (scall == MAT_INITIAL_MATRIX) { 1551 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1552 } 1553 1554 for ( i=0; i<n; i++ ) { 1555 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1556 } 1557 PetscFunctionReturn(0); 1558 } 1559 1560 #undef __FUNC__ 1561 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1562 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1563 { 1564 PetscFunctionBegin; 1565 *bs = 1; 1566 PetscFunctionReturn(0); 1567 } 1568 1569 #undef __FUNC__ 1570 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1571 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1572 { 1573 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1574 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1575 int start, end, *ai, *aj; 1576 BT table; 1577 1578 PetscFunctionBegin; 1579 shift = a->indexshift; 1580 m = a->m; 1581 ai = a->i; 1582 aj = a->j+shift; 1583 1584 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1585 1586 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1587 ierr = BTCreate(m,table); CHKERRQ(ierr); 1588 1589 for ( i=0; i<is_max; i++ ) { 1590 /* Initialize the two local arrays */ 1591 isz = 0; 1592 BTMemzero(m,table); 1593 1594 /* Extract the indices, assume there can be duplicate entries */ 1595 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1596 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1597 1598 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1599 for ( j=0; j<n ; ++j){ 1600 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1601 } 1602 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1603 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1604 1605 k = 0; 1606 for ( j=0; j<ov; j++){ /* for each overlap */ 1607 n = isz; 1608 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1609 row = nidx[k]; 1610 start = ai[row]; 1611 end = ai[row+1]; 1612 for ( l = start; l<end ; l++){ 1613 val = aj[l] + shift; 1614 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1615 } 1616 } 1617 } 1618 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1619 } 1620 BTDestroy(table); 1621 PetscFree(nidx); 1622 PetscFunctionReturn(0); 1623 } 1624 1625 /* -------------------------------------------------------------- */ 1626 #undef __FUNC__ 1627 #define __FUNC__ "MatPermute_SeqAIJ" 1628 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1629 { 1630 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1631 Scalar *vwork; 1632 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1633 int *row,*col,*cnew,j,*lens; 1634 IS icolp,irowp; 1635 1636 PetscFunctionBegin; 1637 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1638 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1639 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1640 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1641 1642 /* determine lengths of permuted rows */ 1643 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1644 for (i=0; i<m; i++ ) { 1645 lens[row[i]] = a->i[i+1] - a->i[i]; 1646 } 1647 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1648 PetscFree(lens); 1649 1650 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1651 for (i=0; i<m; i++) { 1652 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1653 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1654 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1655 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1656 } 1657 PetscFree(cnew); 1658 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1659 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1660 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1661 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1662 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1663 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1664 PetscFunctionReturn(0); 1665 } 1666 1667 #undef __FUNC__ 1668 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1669 int MatPrintHelp_SeqAIJ(Mat A) 1670 { 1671 static int called = 0; 1672 MPI_Comm comm = A->comm; 1673 1674 PetscFunctionBegin; 1675 if (called) {PetscFunctionReturn(0);} else called = 1; 1676 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1677 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1678 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1679 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1680 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1681 #if defined(HAVE_ESSL) 1682 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1683 #endif 1684 PetscFunctionReturn(0); 1685 } 1686 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1687 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1688 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1689 1690 /* -------------------------------------------------------------------*/ 1691 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1692 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1693 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1694 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1695 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1696 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1697 MatLUFactor_SeqAIJ,0, 1698 MatRelax_SeqAIJ, 1699 MatTranspose_SeqAIJ, 1700 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1701 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1702 0,MatAssemblyEnd_SeqAIJ, 1703 MatCompress_SeqAIJ, 1704 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1705 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1706 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1707 MatILUFactorSymbolic_SeqAIJ,0, 1708 0,0, 1709 MatConvertSameType_SeqAIJ,0,0, 1710 MatILUFactor_SeqAIJ,0,0, 1711 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1712 MatGetValues_SeqAIJ,0, 1713 MatPrintHelp_SeqAIJ, 1714 MatScale_SeqAIJ,0,0, 1715 MatILUDTFactor_SeqAIJ, 1716 MatGetBlockSize_SeqAIJ, 1717 MatGetRowIJ_SeqAIJ, 1718 MatRestoreRowIJ_SeqAIJ, 1719 MatGetColumnIJ_SeqAIJ, 1720 MatRestoreColumnIJ_SeqAIJ, 1721 MatFDColoringCreate_SeqAIJ, 1722 MatColoringPatch_SeqAIJ, 1723 0, 1724 MatPermute_SeqAIJ, 1725 0, 1726 0, 1727 0, 1728 0, 1729 MatGetMaps_Petsc}; 1730 1731 extern int MatUseSuperLU_SeqAIJ(Mat); 1732 extern int MatUseEssl_SeqAIJ(Mat); 1733 extern int MatUseDXML_SeqAIJ(Mat); 1734 1735 #undef __FUNC__ 1736 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1737 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1738 { 1739 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1740 int i,nz,n; 1741 1742 PetscFunctionBegin; 1743 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1744 1745 nz = aij->maxnz; 1746 n = aij->n; 1747 for (i=0; i<nz; i++) { 1748 aij->j[i] = indices[i]; 1749 } 1750 aij->nz = nz; 1751 for ( i=0; i<n; i++ ) { 1752 aij->ilen[i] = aij->imax[i]; 1753 } 1754 1755 PetscFunctionReturn(0); 1756 } 1757 1758 #undef __FUNC__ 1759 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1760 /*@ 1761 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1762 in the matrix. 1763 1764 Input Parameters: 1765 + mat - the SeqAIJ matrix 1766 - indices - the column indices 1767 1768 Notes: 1769 This can be called if you have precomputed the nonzero structure of the 1770 matrix and want to provide it to the matrix object to improve the performance 1771 of the MatSetValues() operation. 1772 1773 You MUST have set the correct numbers of nonzeros per row in the call to 1774 MatCreateSeqAIJ(). 1775 1776 MUST be called before any calls to MatSetValues(); 1777 1778 @*/ 1779 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1780 { 1781 int ierr,(*f)(Mat,int *); 1782 1783 PetscFunctionBegin; 1784 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1785 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1786 CHKERRQ(ierr); 1787 if (f) { 1788 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1789 } else { 1790 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1791 } 1792 PetscFunctionReturn(0); 1793 } 1794 1795 #undef __FUNC__ 1796 #define __FUNC__ "MatCreateSeqAIJ" 1797 /*@C 1798 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1799 (the default parallel PETSc format). For good matrix assembly performance 1800 the user should preallocate the matrix storage by setting the parameter nz 1801 (or the array nzz). By setting these parameters accurately, performance 1802 during matrix assembly can be increased by more than a factor of 50. 1803 1804 Collective on MPI_Comm 1805 1806 Input Parameters: 1807 + comm - MPI communicator, set to PETSC_COMM_SELF 1808 . m - number of rows 1809 . n - number of columns 1810 . nz - number of nonzeros per row (same for all rows) 1811 - nzz - array containing the number of nonzeros in the various rows 1812 (possibly different for each row) or PETSC_NULL 1813 1814 Output Parameter: 1815 . A - the matrix 1816 1817 Notes: 1818 The AIJ format (also called the Yale sparse matrix format or 1819 compressed row storage), is fully compatible with standard Fortran 77 1820 storage. That is, the stored row and column indices can begin at 1821 either one (as in Fortran) or zero. See the users' manual for details. 1822 1823 Specify the preallocated storage with either nz or nnz (not both). 1824 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1825 allocation. For large problems you MUST preallocate memory or you 1826 will get TERRIBLE performance, see the users' manual chapter on matrices. 1827 1828 By default, this format uses inodes (identical nodes) when possible, to 1829 improve numerical efficiency of matrix-vector products and solves. We 1830 search for consecutive rows with the same nonzero structure, thereby 1831 reusing matrix information to achieve increased efficiency. 1832 1833 Options Database Keys: 1834 + -mat_aij_no_inode - Do not use inodes 1835 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 1836 - -mat_aij_oneindex - Internally use indexing starting at 1 1837 rather than 0. Note that when calling MatSetValues(), 1838 the user still MUST index entries starting at 0! 1839 1840 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 1841 @*/ 1842 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1843 { 1844 Mat B; 1845 Mat_SeqAIJ *b; 1846 int i, len, ierr, flg,size; 1847 1848 PetscFunctionBegin; 1849 MPI_Comm_size(comm,&size); 1850 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 1851 1852 *A = 0; 1853 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1854 PLogObjectCreate(B); 1855 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1856 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1857 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1858 B->ops->destroy = MatDestroy_SeqAIJ; 1859 B->ops->view = MatView_SeqAIJ; 1860 B->factor = 0; 1861 B->lupivotthreshold = 1.0; 1862 B->mapping = 0; 1863 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 1864 b->ilu_preserve_row_sums = PETSC_FALSE; 1865 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 1866 b->row = 0; 1867 b->col = 0; 1868 b->icol = 0; 1869 b->indexshift = 0; 1870 b->reallocs = 0; 1871 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1872 if (flg) b->indexshift = -1; 1873 1874 b->m = m; B->m = m; B->M = m; 1875 b->n = n; B->n = n; B->N = n; 1876 1877 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1878 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1879 1880 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1881 if (nnz == PETSC_NULL) { 1882 if (nz == PETSC_DEFAULT) nz = 10; 1883 else if (nz <= 0) nz = 1; 1884 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1885 nz = nz*m; 1886 } else { 1887 nz = 0; 1888 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1889 } 1890 1891 /* allocate the matrix space */ 1892 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1893 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1894 b->j = (int *) (b->a + nz); 1895 PetscMemzero(b->j,nz*sizeof(int)); 1896 b->i = b->j + nz; 1897 b->singlemalloc = 1; 1898 1899 b->i[0] = -b->indexshift; 1900 for (i=1; i<m+1; i++) { 1901 b->i[i] = b->i[i-1] + b->imax[i-1]; 1902 } 1903 1904 /* b->ilen will count nonzeros in each row so far. */ 1905 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1906 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1907 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1908 1909 b->nz = 0; 1910 b->maxnz = nz; 1911 b->sorted = 0; 1912 b->roworiented = 1; 1913 b->nonew = 0; 1914 b->diag = 0; 1915 b->solve_work = 0; 1916 b->spptr = 0; 1917 b->inode.node_count = 0; 1918 b->inode.size = 0; 1919 b->inode.limit = 5; 1920 b->inode.max_limit = 5; 1921 B->info.nz_unneeded = (double)b->maxnz; 1922 1923 *A = B; 1924 1925 /* SuperLU is not currently supported through PETSc */ 1926 #if defined(HAVE_SUPERLU) 1927 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1928 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1929 #endif 1930 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1931 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1932 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1933 if (flg) { 1934 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1935 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1936 } 1937 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1938 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1939 1940 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 1941 "MatSeqAIJSetColumnIndices_SeqAIJ", 1942 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 1943 PetscFunctionReturn(0); 1944 } 1945 1946 #undef __FUNC__ 1947 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1948 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1949 { 1950 Mat C; 1951 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1952 int i,len, m = a->m,shift = a->indexshift,ierr; 1953 1954 PetscFunctionBegin; 1955 *B = 0; 1956 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1957 PLogObjectCreate(C); 1958 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1959 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1960 C->ops->destroy = MatDestroy_SeqAIJ; 1961 C->ops->view = MatView_SeqAIJ; 1962 C->factor = A->factor; 1963 c->row = 0; 1964 c->col = 0; 1965 c->icol = 0; 1966 c->indexshift = shift; 1967 C->assembled = PETSC_TRUE; 1968 1969 c->m = C->m = a->m; 1970 c->n = C->n = a->n; 1971 C->M = a->m; 1972 C->N = a->n; 1973 1974 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1975 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1976 for ( i=0; i<m; i++ ) { 1977 c->imax[i] = a->imax[i]; 1978 c->ilen[i] = a->ilen[i]; 1979 } 1980 1981 /* allocate the matrix space */ 1982 c->singlemalloc = 1; 1983 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1984 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1985 c->j = (int *) (c->a + a->i[m] + shift); 1986 c->i = c->j + a->i[m] + shift; 1987 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1988 if (m > 0) { 1989 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1990 if (cpvalues == COPY_VALUES) { 1991 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1992 } 1993 } 1994 1995 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1996 c->sorted = a->sorted; 1997 c->roworiented = a->roworiented; 1998 c->nonew = a->nonew; 1999 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2000 2001 if (a->diag) { 2002 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2003 PLogObjectMemory(C,(m+1)*sizeof(int)); 2004 for ( i=0; i<m; i++ ) { 2005 c->diag[i] = a->diag[i]; 2006 } 2007 } else c->diag = 0; 2008 c->inode.limit = a->inode.limit; 2009 c->inode.max_limit = a->inode.max_limit; 2010 if (a->inode.size){ 2011 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2012 c->inode.node_count = a->inode.node_count; 2013 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2014 } else { 2015 c->inode.size = 0; 2016 c->inode.node_count = 0; 2017 } 2018 c->nz = a->nz; 2019 c->maxnz = a->maxnz; 2020 c->solve_work = 0; 2021 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2022 2023 *B = C; 2024 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2025 "MatSeqAIJSetColumnIndices_SeqAIJ", 2026 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2027 PetscFunctionReturn(0); 2028 } 2029 2030 #undef __FUNC__ 2031 #define __FUNC__ "MatLoad_SeqAIJ" 2032 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2033 { 2034 Mat_SeqAIJ *a; 2035 Mat B; 2036 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2037 MPI_Comm comm; 2038 2039 PetscFunctionBegin; 2040 PetscObjectGetComm((PetscObject) viewer,&comm); 2041 MPI_Comm_size(comm,&size); 2042 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2043 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2044 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2045 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2046 M = header[1]; N = header[2]; nz = header[3]; 2047 2048 if (nz < 0) { 2049 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2050 } 2051 2052 /* read in row lengths */ 2053 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2054 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2055 2056 /* create our matrix */ 2057 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2058 B = *A; 2059 a = (Mat_SeqAIJ *) B->data; 2060 shift = a->indexshift; 2061 2062 /* read in column indices and adjust for Fortran indexing*/ 2063 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2064 if (shift) { 2065 for ( i=0; i<nz; i++ ) { 2066 a->j[i] += 1; 2067 } 2068 } 2069 2070 /* read in nonzero values */ 2071 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2072 2073 /* set matrix "i" values */ 2074 a->i[0] = -shift; 2075 for ( i=1; i<= M; i++ ) { 2076 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2077 a->ilen[i-1] = rowlengths[i-1]; 2078 } 2079 PetscFree(rowlengths); 2080 2081 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2082 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2083 PetscFunctionReturn(0); 2084 } 2085 2086 #undef __FUNC__ 2087 #define __FUNC__ "MatEqual_SeqAIJ" 2088 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2089 { 2090 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2091 2092 PetscFunctionBegin; 2093 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2094 2095 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2096 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2097 (a->indexshift != b->indexshift)) { 2098 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2099 } 2100 2101 /* if the a->i are the same */ 2102 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2103 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2104 } 2105 2106 /* if a->j are the same */ 2107 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2108 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2109 } 2110 2111 /* if a->a are the same */ 2112 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2113 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2114 } 2115 *flg = PETSC_TRUE; 2116 PetscFunctionReturn(0); 2117 2118 } 2119