1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.327 1999/09/27 21:29:41 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 "sys.h" 11 #include "src/mat/impls/aij/seq/aij.h" 12 #include "src/vec/vecimpl.h" 13 #include "src/inline/spops.h" 14 #include "src/inline/dot.h" 15 #include "bitarray.h" 16 17 /* 18 Basic AIJ format ILU based on drop tolerance 19 */ 20 #undef __FUNC__ 21 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 22 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 23 { 24 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 25 26 PetscFunctionBegin; 27 SETERRQ(1,0,"Not implemented"); 28 #if !defined(PETSC_USE_DEBUG) 29 PetscFunctionReturn(0); 30 #endif 31 } 32 33 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 34 35 #undef __FUNC__ 36 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 37 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 38 PetscTruth *done) 39 { 40 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 41 int ierr,i,ishift; 42 43 PetscFunctionBegin; 44 *m = A->m; 45 if (!ia) PetscFunctionReturn(0); 46 ishift = a->indexshift; 47 if (symmetric) { 48 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 49 } else if (oshift == 0 && ishift == -1) { 50 int nz = a->i[a->m]; 51 /* malloc space and subtract 1 from i and j indices */ 52 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) );CHKPTRQ(*ia); 53 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) );CHKPTRQ(*ja); 54 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 55 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 56 } else if (oshift == 1 && ishift == 0) { 57 int nz = a->i[a->m] + 1; 58 /* malloc space and add 1 to i and j indices */ 59 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) );CHKPTRQ(*ia); 60 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) );CHKPTRQ(*ja); 61 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 62 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 63 } else { 64 *ia = a->i; *ja = a->j; 65 } 66 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNC__ 71 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 72 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 73 PetscTruth *done) 74 { 75 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 76 int ishift = a->indexshift,ierr; 77 78 PetscFunctionBegin; 79 if (!ia) PetscFunctionReturn(0); 80 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 81 ierr = PetscFree(*ia);CHKERRQ(ierr); 82 ierr = PetscFree(*ja);CHKERRQ(ierr); 83 } 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNC__ 88 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 89 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 90 PetscTruth *done) 91 { 92 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 93 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 94 int nz = a->i[m]+ishift,row,*jj,mr,col; 95 96 PetscFunctionBegin; 97 *nn = A->n; 98 if (!ia) PetscFunctionReturn(0); 99 if (symmetric) { 100 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 101 } else { 102 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(collengths); 103 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 104 cia = (int *) PetscMalloc( (n+1)*sizeof(int) );CHKPTRQ(cia); 105 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) );CHKPTRQ(cja); 106 jj = a->j; 107 for ( i=0; i<nz; i++ ) { 108 collengths[jj[i] + ishift]++; 109 } 110 cia[0] = oshift; 111 for ( i=0; i<n; i++) { 112 cia[i+1] = cia[i] + collengths[i]; 113 } 114 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 115 jj = a->j; 116 for ( row=0; row<m; row++ ) { 117 mr = a->i[row+1] - a->i[row]; 118 for ( i=0; i<mr; i++ ) { 119 col = *jj++ + ishift; 120 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 121 } 122 } 123 ierr = PetscFree(collengths);CHKERRQ(ierr); 124 *ia = cia; *ja = cja; 125 } 126 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNC__ 131 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 132 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 133 int **ja,PetscTruth *done) 134 { 135 int ierr; 136 137 PetscFunctionBegin; 138 if (!ia) PetscFunctionReturn(0); 139 140 ierr = PetscFree(*ia);CHKERRQ(ierr); 141 ierr = PetscFree(*ja);CHKERRQ(ierr); 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,ierr; 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 (row < 0) continue; 162 #if defined(PETSC_USE_BOPT_g) 163 if (row >= a->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large: row %d max %d",row,a->m); 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 (in[l] < 0) continue; 170 #if defined(PETSC_USE_BOPT_g) 171 if (in[l] >= a->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large: col %d max %d",in[l],a->n); 172 #endif 173 col = in[l] - shift; 174 if (roworiented) { 175 value = v[l + k*n]; 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 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int));CHKERRQ(ierr); 212 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 213 ierr = PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow,len*sizeof(int));CHKERRQ(ierr); 214 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar));CHKERRQ(ierr); 215 ierr = PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow,len*sizeof(Scalar));CHKERRQ(ierr); 216 /* free up old matrix storage */ 217 ierr = PetscFree(a->a);CHKERRQ(ierr); 218 if (!a->singlemalloc) { 219 ierr = PetscFree(a->i);CHKERRQ(ierr); 220 ierr = PetscFree(a->j);CHKERRQ(ierr); 221 } 222 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 223 a->singlemalloc = 1; 224 225 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 226 rmax = imax[row] = imax[row] + CHUNKSIZE; 227 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 228 a->maxnz += CHUNKSIZE; 229 a->reallocs++; 230 } 231 N = nrow++ - 1; a->nz++; 232 /* shift up all the later entries in this row */ 233 for ( ii=N; ii>=i; ii-- ) { 234 rp[ii+1] = rp[ii]; 235 ap[ii+1] = ap[ii]; 236 } 237 rp[i] = col; 238 ap[i] = value; 239 noinsert:; 240 low = i + 1; 241 } 242 ailen[row] = nrow; 243 } 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNC__ 248 #define __FUNC__ "MatGetValues_SeqAIJ" 249 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 250 { 251 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 252 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 253 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 254 Scalar *ap, *aa = a->a, zero = 0.0; 255 256 PetscFunctionBegin; 257 for ( k=0; k<m; k++ ) { /* loop over rows */ 258 row = im[k]; 259 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 260 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 261 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 262 nrow = ailen[row]; 263 for ( l=0; l<n; l++ ) { /* loop over columns */ 264 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 265 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 266 col = in[l] - shift; 267 high = nrow; low = 0; /* assume unsorted */ 268 while (high-low > 5) { 269 t = (low+high)/2; 270 if (rp[t] > col) high = t; 271 else low = t; 272 } 273 for ( i=low; i<high; i++ ) { 274 if (rp[i] > col) break; 275 if (rp[i] == col) { 276 *v++ = ap[i]; 277 goto finished; 278 } 279 } 280 *v++ = zero; 281 finished:; 282 } 283 } 284 PetscFunctionReturn(0); 285 } 286 287 288 #undef __FUNC__ 289 #define __FUNC__ "MatView_SeqAIJ_Binary" 290 int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 291 { 292 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 293 int i, fd, *col_lens, ierr; 294 295 PetscFunctionBegin; 296 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 297 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) );CHKPTRQ(col_lens); 298 col_lens[0] = MAT_COOKIE; 299 col_lens[1] = a->m; 300 col_lens[2] = a->n; 301 col_lens[3] = a->nz; 302 303 /* store lengths of each row and write (including header) to file */ 304 for ( i=0; i<a->m; i++ ) { 305 col_lens[4+i] = a->i[i+1] - a->i[i]; 306 } 307 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1);CHKERRQ(ierr); 308 ierr = PetscFree(col_lens);CHKERRQ(ierr); 309 310 /* store column indices (zero start index) */ 311 if (a->indexshift) { 312 for ( i=0; i<a->nz; i++ ) a->j[i]--; 313 } 314 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0);CHKERRQ(ierr); 315 if (a->indexshift) { 316 for ( i=0; i<a->nz; i++ ) a->j[i]++; 317 } 318 319 /* store nonzero values */ 320 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 321 PetscFunctionReturn(0); 322 } 323 324 #undef __FUNC__ 325 #define __FUNC__ "MatView_SeqAIJ_ASCII" 326 int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 327 { 328 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 329 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 330 FILE *fd; 331 char *outputname; 332 333 PetscFunctionBegin; 334 ierr = ViewerASCIIGetPointer(viewer,&fd);CHKERRQ(ierr); 335 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 336 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 337 if (format == VIEWER_FORMAT_ASCII_INFO) { 338 PetscFunctionReturn(0); 339 } else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 340 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1);CHKERRQ(ierr); 341 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2);CHKERRQ(ierr); 342 if (flg1 || flg2) {ierr = ViewerASCIIPrintf(viewer," not using I-node routines\n");CHKERRQ(ierr);} 343 else {ierr = ViewerASCIIPrintf(viewer," using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr);} 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(PETSC_USE_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 if (outputname) fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 367 else fprintf(fd,"];\n M = spconvert(zzz);\n"); 368 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 369 for ( i=0; i<m; i++ ) { 370 fprintf(fd,"row %d:",i); 371 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 372 #if defined(PETSC_USE_COMPLEX) 373 if (PetscImaginary(a->a[j]) > 0.0 && PetscReal(a->a[j]) != 0.0) 374 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 375 else if (PetscImaginary(a->a[j]) < 0.0 && PetscReal(a->a[j]) != 0.0) 376 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 377 else if (PetscReal(a->a[j]) != 0.0) 378 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(a->a[j])); 379 #else 380 if (a->a[j] != 0.0) fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 381 #endif 382 } 383 fprintf(fd,"\n"); 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(PETSC_USE_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 ierr = PetscFree(sptr);CHKERRQ(ierr); 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(PETSC_USE_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 if (format == VIEWER_FORMAT_ASCII_DENSE) { 433 int cnt = 0,jcnt; 434 Scalar value; 435 436 for ( i=0; i<m; i++ ) { 437 jcnt = 0; 438 for ( j=0; j<a->n; j++ ) { 439 if ( jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 440 value = a->a[cnt++]; 441 jcnt++; 442 } else { 443 value = 0.0; 444 } 445 #if defined(PETSC_USE_COMPLEX) 446 fprintf(fd," %7.5e+%7.5e i ",PetscReal(value),PetscImaginary(value)); 447 #else 448 fprintf(fd," %7.5e ",value); 449 #endif 450 } 451 fprintf(fd,"\n"); 452 } 453 } else { 454 for ( i=0; i<m; i++ ) { 455 fprintf(fd,"row %d:",i); 456 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 457 #if defined(PETSC_USE_COMPLEX) 458 if (PetscImaginary(a->a[j]) > 0.0) { 459 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 460 } else if (PetscImaginary(a->a[j]) < 0.0) { 461 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 462 } else { 463 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(a->a[j])); 464 } 465 #else 466 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 467 #endif 468 } 469 fprintf(fd,"\n"); 470 } 471 } 472 fflush(fd); 473 PetscFunctionReturn(0); 474 } 475 476 #undef __FUNC__ 477 #define __FUNC__ "MatView_SeqAIJ_Draw_Zoom" 478 int MatView_SeqAIJ_Draw_Zoom(Draw draw,void *Aa) 479 { 480 Mat A = (Mat) Aa; 481 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 482 int ierr, i,j, m = a->m, shift = a->indexshift,color,rank; 483 int format; 484 double xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 485 Viewer viewer; 486 MPI_Comm comm; 487 488 PetscFunctionBegin; 489 /* 490 This is nasty. If this is called from an originally parallel matrix 491 then all processes call this, but only the first has the matrix so the 492 rest should return immediately. 493 */ 494 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 495 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 496 if (rank) PetscFunctionReturn(0); 497 498 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*) &viewer);CHKERRQ(ierr); 499 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 500 501 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 502 /* loop over matrix elements drawing boxes */ 503 504 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 505 /* Blue for negative, Cyan for zero and Red for positive */ 506 color = DRAW_BLUE; 507 for ( i=0; i<m; i++ ) { 508 y_l = m - i - 1.0; y_r = y_l + 1.0; 509 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 510 x_l = a->j[j] + shift; x_r = x_l + 1.0; 511 #if defined(PETSC_USE_COMPLEX) 512 if (PetscReal(a->a[j]) >= 0.) continue; 513 #else 514 if (a->a[j] >= 0.) continue; 515 #endif 516 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 517 } 518 } 519 color = DRAW_CYAN; 520 for ( i=0; i<m; i++ ) { 521 y_l = m - i - 1.0; y_r = y_l + 1.0; 522 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 523 x_l = a->j[j] + shift; x_r = x_l + 1.0; 524 if (a->a[j] != 0.) continue; 525 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 526 } 527 } 528 color = DRAW_RED; 529 for ( i=0; i<m; i++ ) { 530 y_l = m - i - 1.0; y_r = y_l + 1.0; 531 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 532 x_l = a->j[j] + shift; x_r = x_l + 1.0; 533 #if defined(PETSC_USE_COMPLEX) 534 if (PetscReal(a->a[j]) <= 0.) continue; 535 #else 536 if (a->a[j] <= 0.) continue; 537 #endif 538 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 539 } 540 } 541 } else { 542 /* use contour shading to indicate magnitude of values */ 543 /* first determine max of all nonzero values */ 544 int nz = a->nz,count; 545 Draw popup; 546 double scale; 547 548 for ( i=0; i<nz; i++ ) { 549 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 550 } 551 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 552 ierr = DrawGetPopup(draw,&popup);CHKERRQ(ierr); 553 ierr = DrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr); 554 count = 0; 555 for ( i=0; i<m; i++ ) { 556 y_l = m - i - 1.0; y_r = y_l + 1.0; 557 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 558 x_l = a->j[j] + shift; x_r = x_l + 1.0; 559 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 560 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 561 count++; 562 } 563 } 564 } 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNC__ 569 #define __FUNC__ "MatView_SeqAIJ_Draw" 570 int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 571 { 572 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 573 int ierr; 574 Draw draw; 575 double xr,yr,xl,yl,h,w; 576 PetscTruth isnull; 577 578 PetscFunctionBegin; 579 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 580 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); 581 if (isnull) PetscFunctionReturn(0); 582 583 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 584 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 585 xr += w; yr += h; xl = -w; yl = -h; 586 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 587 ierr = DrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 588 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 589 PetscFunctionReturn(0); 590 } 591 592 #undef __FUNC__ 593 #define __FUNC__ "MatView_SeqAIJ" 594 int MatView_SeqAIJ(Mat A,Viewer viewer) 595 { 596 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 597 int ierr; 598 599 PetscFunctionBegin; 600 if (PetscTypeCompare(viewer,SOCKET_VIEWER)) { 601 ierr = ViewerSocketPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 602 } else if (PetscTypeCompare(viewer,ASCII_VIEWER)){ 603 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 604 } else if (PetscTypeCompare(viewer,BINARY_VIEWER)) { 605 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 606 } else if (PetscTypeCompare(viewer,DRAW_VIEWER)) { 607 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 608 } else { 609 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 610 } 611 PetscFunctionReturn(0); 612 } 613 614 extern int Mat_AIJ_CheckInode(Mat); 615 #undef __FUNC__ 616 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 617 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 618 { 619 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 620 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 621 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax = 0; 622 Scalar *aa = a->a, *ap; 623 624 PetscFunctionBegin; 625 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 626 627 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 628 for ( i=1; i<m; i++ ) { 629 /* move each row back by the amount of empty slots (fshift) before it*/ 630 fshift += imax[i-1] - ailen[i-1]; 631 rmax = PetscMax(rmax,ailen[i]); 632 if (fshift) { 633 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 634 N = ailen[i]; 635 for ( j=0; j<N; j++ ) { 636 ip[j-fshift] = ip[j]; 637 ap[j-fshift] = ap[j]; 638 } 639 } 640 ai[i] = ai[i-1] + ailen[i-1]; 641 } 642 if (m) { 643 fshift += imax[m-1] - ailen[m-1]; 644 ai[m] = ai[m-1] + ailen[m-1]; 645 } 646 /* reset ilen and imax for each row */ 647 for ( i=0; i<m; i++ ) { 648 ailen[i] = imax[i] = ai[i+1] - ai[i]; 649 } 650 a->nz = ai[m] + shift; 651 652 /* diagonals may have moved, so kill the diagonal pointers */ 653 if (fshift && a->diag) { 654 ierr = PetscFree(a->diag);CHKERRQ(ierr); 655 PLogObjectMemory(A,-(m+1)*sizeof(int)); 656 a->diag = 0; 657 } 658 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 659 m,a->n,fshift,a->nz); 660 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n", 661 a->reallocs); 662 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 663 a->reallocs = 0; 664 A->info.nz_unneeded = (double)fshift; 665 666 /* check out for identical nodes. If found, use inode functions */ 667 ierr = Mat_AIJ_CheckInode(A);CHKERRQ(ierr); 668 PetscFunctionReturn(0); 669 } 670 671 #undef __FUNC__ 672 #define __FUNC__ "MatZeroEntries_SeqAIJ" 673 int MatZeroEntries_SeqAIJ(Mat A) 674 { 675 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 676 int ierr; 677 678 PetscFunctionBegin; 679 ierr = PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 680 PetscFunctionReturn(0); 681 } 682 683 #undef __FUNC__ 684 #define __FUNC__ "MatDestroy_SeqAIJ" 685 int MatDestroy_SeqAIJ(Mat A) 686 { 687 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 688 int ierr; 689 690 PetscFunctionBegin; 691 692 if (A->mapping) { 693 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 694 } 695 if (A->bmapping) { 696 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 697 } 698 if (A->rmap) { 699 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 700 } 701 if (A->cmap) { 702 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 703 } 704 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 705 #if defined(PETSC_USE_LOG) 706 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 707 #endif 708 ierr = PetscFree(a->a);CHKERRQ(ierr); 709 if (!a->singlemalloc) { 710 ierr = PetscFree(a->i);CHKERRQ(ierr); 711 ierr = PetscFree(a->j);CHKERRQ(ierr); 712 } 713 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 714 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 715 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 716 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 717 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 718 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 719 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 720 ierr = PetscFree(a);CHKERRQ(ierr); 721 722 PLogObjectDestroy(A); 723 PetscHeaderDestroy(A); 724 PetscFunctionReturn(0); 725 } 726 727 #undef __FUNC__ 728 #define __FUNC__ "MatCompress_SeqAIJ" 729 int MatCompress_SeqAIJ(Mat A) 730 { 731 PetscFunctionBegin; 732 PetscFunctionReturn(0); 733 } 734 735 #undef __FUNC__ 736 #define __FUNC__ "MatSetOption_SeqAIJ" 737 int MatSetOption_SeqAIJ(Mat A,MatOption op) 738 { 739 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 740 741 PetscFunctionBegin; 742 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 743 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 744 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 745 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 746 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 747 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 748 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 749 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 750 else if (op == MAT_ROWS_SORTED || 751 op == MAT_ROWS_UNSORTED || 752 op == MAT_SYMMETRIC || 753 op == MAT_STRUCTURALLY_SYMMETRIC || 754 op == MAT_YES_NEW_DIAGONALS || 755 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 756 op == MAT_USE_HASH_TABLE) 757 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 758 else if (op == MAT_NO_NEW_DIAGONALS) { 759 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 760 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 761 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 762 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 763 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 764 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 765 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 766 PetscFunctionReturn(0); 767 } 768 769 #undef __FUNC__ 770 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 771 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 772 { 773 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 774 int i,j, n,shift = a->indexshift,ierr; 775 Scalar *x, zero = 0.0; 776 777 PetscFunctionBegin; 778 ierr = VecSet(&zero,v);CHKERRQ(ierr); 779 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 780 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 781 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 782 for ( i=0; i<a->m; i++ ) { 783 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 784 if (a->j[j]+shift == i) { 785 x[i] = a->a[j]; 786 break; 787 } 788 } 789 } 790 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 791 PetscFunctionReturn(0); 792 } 793 794 /* -------------------------------------------------------*/ 795 /* Should check that shapes of vectors and matrices match */ 796 /* -------------------------------------------------------*/ 797 #undef __FUNC__ 798 #define __FUNC__ "MatMultTrans_SeqAIJ" 799 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 800 { 801 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 802 Scalar *x, *y, *v, alpha; 803 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 804 805 PetscFunctionBegin; 806 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 807 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 808 ierr = PetscMemzero(y,a->n*sizeof(Scalar));CHKERRQ(ierr); 809 y = y + shift; /* shift for Fortran start by 1 indexing */ 810 for ( i=0; i<m; i++ ) { 811 idx = a->j + a->i[i] + shift; 812 v = a->a + a->i[i] + shift; 813 n = a->i[i+1] - a->i[i]; 814 alpha = x[i]; 815 while (n-->0) {y[*idx++] += alpha * *v++;} 816 } 817 PLogFlops(2*a->nz - a->n); 818 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 819 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 820 PetscFunctionReturn(0); 821 } 822 823 #undef __FUNC__ 824 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 825 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 826 { 827 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 828 Scalar *x, *y, *v, alpha; 829 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 830 831 PetscFunctionBegin; 832 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 833 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 834 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 835 y = y + shift; /* shift for Fortran start by 1 indexing */ 836 for ( i=0; i<m; i++ ) { 837 idx = a->j + a->i[i] + shift; 838 v = a->a + a->i[i] + shift; 839 n = a->i[i+1] - a->i[i]; 840 alpha = x[i]; 841 while (n-->0) {y[*idx++] += alpha * *v++;} 842 } 843 PLogFlops(2*a->nz); 844 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 845 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 846 PetscFunctionReturn(0); 847 } 848 849 #undef __FUNC__ 850 #define __FUNC__ "MatMult_SeqAIJ" 851 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 852 { 853 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 854 Scalar *x, *y, *v, sum; 855 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 856 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 857 int n, i, jrow,j; 858 #endif 859 860 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 861 #pragma disjoint(*x,*y,*v) 862 #endif 863 864 PetscFunctionBegin; 865 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 866 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 867 x = x + shift; /* shift for Fortran start by 1 indexing */ 868 idx = a->j; 869 v = a->a; 870 ii = a->i; 871 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 872 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 873 #else 874 v += shift; /* shift for Fortran start by 1 indexing */ 875 idx += shift; 876 for ( i=0; i<m; i++ ) { 877 jrow = ii[i]; 878 n = ii[i+1] - jrow; 879 sum = 0.0; 880 for ( j=0; j<n; j++) { 881 sum += v[jrow]*x[idx[jrow]]; jrow++; 882 } 883 y[i] = sum; 884 } 885 #endif 886 PLogFlops(2*a->nz - m); 887 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 888 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 889 PetscFunctionReturn(0); 890 } 891 892 #undef __FUNC__ 893 #define __FUNC__ "MatMultAdd_SeqAIJ" 894 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 895 { 896 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 897 Scalar *x, *y, *z, *v, sum; 898 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 899 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 900 int n,i,jrow,j; 901 #endif 902 903 PetscFunctionBegin; 904 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 905 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 906 if (zz != yy) { 907 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 908 } else { 909 z = y; 910 } 911 x = x + shift; /* shift for Fortran start by 1 indexing */ 912 idx = a->j; 913 v = a->a; 914 ii = a->i; 915 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 916 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 917 #else 918 v += shift; /* shift for Fortran start by 1 indexing */ 919 idx += shift; 920 for ( i=0; i<m; i++ ) { 921 jrow = ii[i]; 922 n = ii[i+1] - jrow; 923 sum = y[i]; 924 for ( j=0; j<n; j++) { 925 sum += v[jrow]*x[idx[jrow]]; jrow++; 926 } 927 z[i] = sum; 928 } 929 #endif 930 PLogFlops(2*a->nz); 931 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 932 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 933 if (zz != yy) { 934 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 935 } 936 PetscFunctionReturn(0); 937 } 938 939 /* 940 Adds diagonal pointers to sparse matrix structure. 941 */ 942 #undef __FUNC__ 943 #define __FUNC__ "MatMarkDiag_SeqAIJ" 944 int MatMarkDiag_SeqAIJ(Mat A) 945 { 946 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 947 int i,j, *diag, m = a->m,shift = a->indexshift; 948 949 PetscFunctionBegin; 950 diag = (int *) PetscMalloc( (m+1)*sizeof(int));CHKPTRQ(diag); 951 PLogObjectMemory(A,(m+1)*sizeof(int)); 952 for ( i=0; i<a->m; i++ ) { 953 diag[i] = a->i[i+1]; 954 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 955 if (a->j[j]+shift == i) { 956 diag[i] = j - shift; 957 break; 958 } 959 } 960 } 961 a->diag = diag; 962 PetscFunctionReturn(0); 963 } 964 965 /* 966 Checks for missing diagonals 967 */ 968 #undef __FUNC__ 969 #define __FUNC__ "MatMissingDiag_SeqAIJ" 970 int MatMissingDiag_SeqAIJ(Mat A) 971 { 972 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 973 int *diag = a->diag, *jj = a->j,i,shift = a->indexshift; 974 975 PetscFunctionBegin; 976 for ( i=0; i<a->m; i++ ) { 977 if (jj[diag[i]+shift] != i-shift) { 978 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 979 } 980 } 981 PetscFunctionReturn(0); 982 } 983 984 #undef __FUNC__ 985 #define __FUNC__ "MatRelax_SeqAIJ" 986 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag,double fshift,int its,Vec xx) 987 { 988 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 989 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t=0,scale,*ts, *xb,*idiag=0; 990 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 991 992 PetscFunctionBegin; 993 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 994 if (xx != bb) { 995 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 996 } else { 997 b = x; 998 } 999 1000 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 1001 diag = a->diag; 1002 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1003 if (flag == SOR_APPLY_UPPER) { 1004 /* apply ( U + D/omega) to the vector */ 1005 bs = b + shift; 1006 for ( i=0; i<m; i++ ) { 1007 d = fshift + a->a[diag[i] + shift]; 1008 n = a->i[i+1] - diag[i] - 1; 1009 PLogFlops(2*n-1); 1010 idx = a->j + diag[i] + (!shift); 1011 v = a->a + diag[i] + (!shift); 1012 sum = b[i]*d/omega; 1013 SPARSEDENSEDOT(sum,bs,v,idx,n); 1014 x[i] = sum; 1015 } 1016 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1017 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1018 PetscFunctionReturn(0); 1019 } 1020 1021 /* setup workspace for Eisenstat */ 1022 if (flag & SOR_EISENSTAT) { 1023 if (!a->idiag) { 1024 a->idiag = (Scalar *) PetscMalloc(2*m*sizeof(Scalar));CHKPTRQ(a->idiag); 1025 a->ssor = a->idiag + m; 1026 v = a->a; 1027 for ( i=0; i<m; i++ ) { a->idiag[i] = 1.0/v[diag[i]];} 1028 } 1029 t = a->ssor; 1030 idiag = a->idiag; 1031 } 1032 /* Let A = L + U + D; where L is lower trianglar, 1033 U is upper triangular, E is diagonal; This routine applies 1034 1035 (L + E)^{-1} A (U + E)^{-1} 1036 1037 to a vector efficiently using Eisenstat's trick. This is for 1038 the case of SSOR preconditioner, so E is D/omega where omega 1039 is the relaxation factor. 1040 */ 1041 1042 if (flag == SOR_APPLY_LOWER) { 1043 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 1044 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1045 /* special case for omega = 1.0 saves flops and some integer ops */ 1046 Scalar *v2; 1047 1048 v2 = a->a; 1049 /* x = (E + U)^{-1} b */ 1050 for ( i=m-1; i>=0; i-- ) { 1051 n = a->i[i+1] - diag[i] - 1; 1052 idx = a->j + diag[i] + 1; 1053 v = a->a + diag[i] + 1; 1054 sum = b[i]; 1055 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1056 x[i] = sum*idiag[i]; 1057 1058 /* t = b - (2*E - D)x */ 1059 t[i] = b[i] - (v2[diag[i]])*x[i]; 1060 } 1061 1062 /* t = (E + L)^{-1}t */ 1063 diag = a->diag; 1064 for ( i=0; i<m; i++ ) { 1065 n = diag[i] - a->i[i]; 1066 idx = a->j + a->i[i]; 1067 v = a->a + a->i[i]; 1068 sum = t[i]; 1069 SPARSEDENSEMDOT(sum,t,v,idx,n); 1070 t[i] = sum*idiag[i]; 1071 1072 /* x = x + t */ 1073 x[i] += t[i]; 1074 } 1075 1076 PLogFlops(3*m-1 + 2*a->nz); 1077 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1078 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1079 PetscFunctionReturn(0); 1080 } else if (flag & SOR_EISENSTAT) { 1081 /* Let A = L + U + D; where L is lower trianglar, 1082 U is upper triangular, E is diagonal; This routine applies 1083 1084 (L + E)^{-1} A (U + E)^{-1} 1085 1086 to a vector efficiently using Eisenstat's trick. This is for 1087 the case of SSOR preconditioner, so E is D/omega where omega 1088 is the relaxation factor. 1089 */ 1090 scale = (2.0/omega) - 1.0; 1091 1092 /* x = (E + U)^{-1} b */ 1093 for ( i=m-1; i>=0; i-- ) { 1094 d = fshift + a->a[diag[i] + shift]; 1095 n = a->i[i+1] - diag[i] - 1; 1096 idx = a->j + diag[i] + (!shift); 1097 v = a->a + diag[i] + (!shift); 1098 sum = b[i]; 1099 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1100 x[i] = omega*(sum/d); 1101 } 1102 1103 /* t = b - (2*E - D)x */ 1104 v = a->a; 1105 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1106 1107 /* t = (E + L)^{-1}t */ 1108 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1109 diag = a->diag; 1110 for ( i=0; i<m; i++ ) { 1111 d = fshift + a->a[diag[i]+shift]; 1112 n = diag[i] - a->i[i]; 1113 idx = a->j + a->i[i] + shift; 1114 v = a->a + a->i[i] + shift; 1115 sum = t[i]; 1116 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1117 t[i] = omega*(sum/d); 1118 /* x = x + t */ 1119 x[i] += t[i]; 1120 } 1121 1122 PLogFlops(6*m-1 + 2*a->nz); 1123 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1124 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1125 PetscFunctionReturn(0); 1126 } 1127 if (flag & SOR_ZERO_INITIAL_GUESS) { 1128 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1129 for ( i=0; i<m; i++ ) { 1130 d = fshift + a->a[diag[i]+shift]; 1131 n = diag[i] - a->i[i]; 1132 PLogFlops(2*n-1); 1133 idx = a->j + a->i[i] + shift; 1134 v = a->a + a->i[i] + shift; 1135 sum = b[i]; 1136 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1137 x[i] = omega*(sum/d); 1138 } 1139 xb = x; 1140 } else xb = b; 1141 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1142 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1143 for ( i=0; i<m; i++ ) { 1144 x[i] *= a->a[diag[i]+shift]; 1145 } 1146 PLogFlops(m); 1147 } 1148 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1149 for ( i=m-1; i>=0; i-- ) { 1150 d = fshift + a->a[diag[i] + shift]; 1151 n = a->i[i+1] - diag[i] - 1; 1152 PLogFlops(2*n-1); 1153 idx = a->j + diag[i] + (!shift); 1154 v = a->a + diag[i] + (!shift); 1155 sum = xb[i]; 1156 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1157 x[i] = omega*(sum/d); 1158 } 1159 } 1160 its--; 1161 } 1162 while (its--) { 1163 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1164 for ( i=0; i<m; i++ ) { 1165 d = fshift + a->a[diag[i]+shift]; 1166 n = a->i[i+1] - a->i[i]; 1167 PLogFlops(2*n-1); 1168 idx = a->j + a->i[i] + shift; 1169 v = a->a + a->i[i] + shift; 1170 sum = b[i]; 1171 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1172 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1173 } 1174 } 1175 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1176 for ( i=m-1; i>=0; i-- ) { 1177 d = fshift + a->a[diag[i] + shift]; 1178 n = a->i[i+1] - a->i[i]; 1179 PLogFlops(2*n-1); 1180 idx = a->j + a->i[i] + shift; 1181 v = a->a + a->i[i] + shift; 1182 sum = b[i]; 1183 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1184 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1185 } 1186 } 1187 } 1188 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1189 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1190 PetscFunctionReturn(0); 1191 } 1192 1193 #undef __FUNC__ 1194 #define __FUNC__ "MatGetInfo_SeqAIJ" 1195 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1196 { 1197 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1198 1199 PetscFunctionBegin; 1200 info->rows_global = (double)a->m; 1201 info->columns_global = (double)a->n; 1202 info->rows_local = (double)a->m; 1203 info->columns_local = (double)a->n; 1204 info->block_size = 1.0; 1205 info->nz_allocated = (double)a->maxnz; 1206 info->nz_used = (double)a->nz; 1207 info->nz_unneeded = (double)(a->maxnz - a->nz); 1208 info->assemblies = (double)A->num_ass; 1209 info->mallocs = (double)a->reallocs; 1210 info->memory = A->mem; 1211 if (A->factor) { 1212 info->fill_ratio_given = A->info.fill_ratio_given; 1213 info->fill_ratio_needed = A->info.fill_ratio_needed; 1214 info->factor_mallocs = A->info.factor_mallocs; 1215 } else { 1216 info->fill_ratio_given = 0; 1217 info->fill_ratio_needed = 0; 1218 info->factor_mallocs = 0; 1219 } 1220 PetscFunctionReturn(0); 1221 } 1222 1223 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1224 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1225 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1226 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1227 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1228 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1229 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1230 1231 #undef __FUNC__ 1232 #define __FUNC__ "MatZeroRows_SeqAIJ" 1233 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1234 { 1235 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1236 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1237 1238 PetscFunctionBegin; 1239 ierr = ISGetSize(is,&N);CHKERRQ(ierr); 1240 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1241 if (diag) { 1242 for ( i=0; i<N; i++ ) { 1243 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1244 if (a->ilen[rows[i]] > 0) { 1245 a->ilen[rows[i]] = 1; 1246 a->a[a->i[rows[i]]+shift] = *diag; 1247 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1248 } else { /* in case row was completely empty */ 1249 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1250 } 1251 } 1252 } else { 1253 for ( i=0; i<N; i++ ) { 1254 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1255 a->ilen[rows[i]] = 0; 1256 } 1257 } 1258 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1259 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262 1263 #undef __FUNC__ 1264 #define __FUNC__ "MatGetSize_SeqAIJ" 1265 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1266 { 1267 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1268 1269 PetscFunctionBegin; 1270 if (m) *m = a->m; 1271 if (n) *n = a->n; 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNC__ 1276 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1277 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1278 { 1279 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1280 1281 PetscFunctionBegin; 1282 *m = 0; *n = a->m; 1283 PetscFunctionReturn(0); 1284 } 1285 1286 #undef __FUNC__ 1287 #define __FUNC__ "MatGetRow_SeqAIJ" 1288 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1289 { 1290 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1291 int *itmp,i,shift = a->indexshift; 1292 1293 PetscFunctionBegin; 1294 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1295 1296 *nz = a->i[row+1] - a->i[row]; 1297 if (v) *v = a->a + a->i[row] + shift; 1298 if (idx) { 1299 itmp = a->j + a->i[row] + shift; 1300 if (*nz && shift) { 1301 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) );CHKPTRQ(*idx); 1302 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1303 } else if (*nz) { 1304 *idx = itmp; 1305 } 1306 else *idx = 0; 1307 } 1308 PetscFunctionReturn(0); 1309 } 1310 1311 #undef __FUNC__ 1312 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1313 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1314 { 1315 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1316 int ierr; 1317 1318 PetscFunctionBegin; 1319 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1320 PetscFunctionReturn(0); 1321 } 1322 1323 #undef __FUNC__ 1324 #define __FUNC__ "MatNorm_SeqAIJ" 1325 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1326 { 1327 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1328 Scalar *v = a->a; 1329 double sum = 0.0; 1330 int i, j,shift = a->indexshift,ierr; 1331 1332 PetscFunctionBegin; 1333 if (type == NORM_FROBENIUS) { 1334 for (i=0; i<a->nz; i++ ) { 1335 #if defined(PETSC_USE_COMPLEX) 1336 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1337 #else 1338 sum += (*v)*(*v); v++; 1339 #endif 1340 } 1341 *norm = sqrt(sum); 1342 } else if (type == NORM_1) { 1343 double *tmp; 1344 int *jj = a->j; 1345 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) );CHKPTRQ(tmp); 1346 ierr = PetscMemzero(tmp,a->n*sizeof(double));CHKERRQ(ierr); 1347 *norm = 0.0; 1348 for ( j=0; j<a->nz; j++ ) { 1349 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1350 } 1351 for ( j=0; j<a->n; j++ ) { 1352 if (tmp[j] > *norm) *norm = tmp[j]; 1353 } 1354 ierr = PetscFree(tmp);CHKERRQ(ierr); 1355 } else if (type == NORM_INFINITY) { 1356 *norm = 0.0; 1357 for ( j=0; j<a->m; j++ ) { 1358 v = a->a + a->i[j] + shift; 1359 sum = 0.0; 1360 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1361 sum += PetscAbsScalar(*v); v++; 1362 } 1363 if (sum > *norm) *norm = sum; 1364 } 1365 } else { 1366 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1367 } 1368 PetscFunctionReturn(0); 1369 } 1370 1371 #undef __FUNC__ 1372 #define __FUNC__ "MatTranspose_SeqAIJ" 1373 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1374 { 1375 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1376 Mat C; 1377 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1378 int shift = a->indexshift; 1379 Scalar *array = a->a; 1380 1381 PetscFunctionBegin; 1382 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1383 col = (int *) PetscMalloc((1+a->n)*sizeof(int));CHKPTRQ(col); 1384 ierr = PetscMemzero(col,(1+a->n)*sizeof(int));CHKERRQ(ierr); 1385 if (shift) { 1386 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1387 } 1388 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1389 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C);CHKERRQ(ierr); 1390 ierr = PetscFree(col);CHKERRQ(ierr); 1391 for ( i=0; i<m; i++ ) { 1392 len = ai[i+1]-ai[i]; 1393 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1394 array += len; aj += len; 1395 } 1396 if (shift) { 1397 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1398 } 1399 1400 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1401 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1402 1403 if (B != PETSC_NULL) { 1404 *B = C; 1405 } else { 1406 PetscOps *Abops; 1407 MatOps Aops; 1408 1409 /* This isn't really an in-place transpose */ 1410 ierr = PetscFree(a->a);CHKERRQ(ierr); 1411 if (!a->singlemalloc) { 1412 ierr = PetscFree(a->i);CHKERRQ(ierr); 1413 ierr = PetscFree(a->j); 1414 } 1415 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 1416 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 1417 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 1418 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 1419 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 1420 ierr = PetscFree(a);CHKERRQ(ierr); 1421 1422 1423 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1424 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1425 1426 /* 1427 This is horrible, horrible code. We need to keep the 1428 the bops and ops Structures, copy everything from C 1429 including the function pointers.. 1430 */ 1431 ierr = PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1432 ierr = PetscMemcpy(A->bops,C->bops,sizeof(PetscOps));CHKERRQ(ierr); 1433 Abops = A->bops; 1434 Aops = A->ops; 1435 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 1436 A->bops = Abops; 1437 A->ops = Aops; 1438 A->qlist = 0; 1439 1440 PetscHeaderDestroy(C); 1441 } 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #undef __FUNC__ 1446 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1447 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1448 { 1449 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1450 Scalar *l,*r,x,*v; 1451 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1452 1453 PetscFunctionBegin; 1454 if (ll) { 1455 /* The local size is used so that VecMPI can be passed to this routine 1456 by MatDiagonalScale_MPIAIJ */ 1457 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1458 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1459 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1460 v = a->a; 1461 for ( i=0; i<m; i++ ) { 1462 x = l[i]; 1463 M = a->i[i+1] - a->i[i]; 1464 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1465 } 1466 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1467 PLogFlops(nz); 1468 } 1469 if (rr) { 1470 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1471 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1472 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1473 v = a->a; jj = a->j; 1474 for ( i=0; i<nz; i++ ) { 1475 (*v++) *= r[*jj++ + shift]; 1476 } 1477 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1478 PLogFlops(nz); 1479 } 1480 PetscFunctionReturn(0); 1481 } 1482 1483 #undef __FUNC__ 1484 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1485 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1486 { 1487 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1488 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1489 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1490 register int sum,lensi; 1491 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1492 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1493 Scalar *a_new,*mat_a; 1494 Mat C; 1495 PetscTruth stride; 1496 1497 PetscFunctionBegin; 1498 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1499 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1500 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1501 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1502 1503 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1504 ierr = ISGetSize(isrow,&nrows);CHKERRQ(ierr); 1505 ierr = ISGetSize(iscol,&ncols);CHKERRQ(ierr); 1506 1507 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1508 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1509 if (stride && step == 1) { 1510 /* special case of contiguous rows */ 1511 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int));CHKPTRQ(lens); 1512 starts = lens + ncols; 1513 /* loop over new rows determining lens and starting points */ 1514 for (i=0; i<nrows; i++) { 1515 kstart = ai[irow[i]]+shift; 1516 kend = kstart + ailen[irow[i]]; 1517 for ( k=kstart; k<kend; k++ ) { 1518 if (aj[k]+shift >= first) { 1519 starts[i] = k; 1520 break; 1521 } 1522 } 1523 sum = 0; 1524 while (k < kend) { 1525 if (aj[k++]+shift >= first+ncols) break; 1526 sum++; 1527 } 1528 lens[i] = sum; 1529 } 1530 /* create submatrix */ 1531 if (scall == MAT_REUSE_MATRIX) { 1532 int n_cols,n_rows; 1533 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1534 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1535 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1536 C = *B; 1537 } else { 1538 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1539 } 1540 c = (Mat_SeqAIJ*) C->data; 1541 1542 /* loop over rows inserting into submatrix */ 1543 a_new = c->a; 1544 j_new = c->j; 1545 i_new = c->i; 1546 i_new[0] = -shift; 1547 for (i=0; i<nrows; i++) { 1548 ii = starts[i]; 1549 lensi = lens[i]; 1550 for ( k=0; k<lensi; k++ ) { 1551 *j_new++ = aj[ii+k] - first; 1552 } 1553 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar));CHKERRQ(ierr); 1554 a_new += lensi; 1555 i_new[i+1] = i_new[i] + lensi; 1556 c->ilen[i] = lensi; 1557 } 1558 ierr = PetscFree(lens);CHKERRQ(ierr); 1559 } else { 1560 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1561 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int));CHKPTRQ(smap); 1562 ssmap = smap + shift; 1563 lens = (int *) PetscMalloc((1+nrows)*sizeof(int));CHKPTRQ(lens); 1564 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1565 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1566 /* determine lens of each row */ 1567 for (i=0; i<nrows; i++) { 1568 kstart = ai[irow[i]]+shift; 1569 kend = kstart + a->ilen[irow[i]]; 1570 lens[i] = 0; 1571 for ( k=kstart; k<kend; k++ ) { 1572 if (ssmap[aj[k]]) { 1573 lens[i]++; 1574 } 1575 } 1576 } 1577 /* Create and fill new matrix */ 1578 if (scall == MAT_REUSE_MATRIX) { 1579 c = (Mat_SeqAIJ *)((*B)->data); 1580 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1581 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1582 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1583 } 1584 ierr = PetscMemzero(c->ilen,c->m*sizeof(int));CHKERRQ(ierr); 1585 C = *B; 1586 } else { 1587 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1588 } 1589 c = (Mat_SeqAIJ *)(C->data); 1590 for (i=0; i<nrows; i++) { 1591 row = irow[i]; 1592 kstart = ai[row]+shift; 1593 kend = kstart + a->ilen[row]; 1594 mat_i = c->i[i]+shift; 1595 mat_j = c->j + mat_i; 1596 mat_a = c->a + mat_i; 1597 mat_ilen = c->ilen + i; 1598 for ( k=kstart; k<kend; k++ ) { 1599 if ((tcol=ssmap[a->j[k]])) { 1600 *mat_j++ = tcol - (!shift); 1601 *mat_a++ = a->a[k]; 1602 (*mat_ilen)++; 1603 1604 } 1605 } 1606 } 1607 /* Free work space */ 1608 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1609 ierr = PetscFree(smap);CHKERRQ(ierr); 1610 ierr = PetscFree(lens);CHKERRQ(ierr); 1611 } 1612 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1613 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1614 1615 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1616 *B = C; 1617 PetscFunctionReturn(0); 1618 } 1619 1620 /* 1621 */ 1622 #undef __FUNC__ 1623 #define __FUNC__ "MatILUFactor_SeqAIJ" 1624 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1625 { 1626 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1627 int ierr; 1628 Mat outA; 1629 PetscTruth row_identity, col_identity; 1630 1631 PetscFunctionBegin; 1632 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,0,"Only levels=0 supported for in-place ilu"); 1633 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1634 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1635 if (!row_identity || !col_identity) { 1636 SETERRQ(1,1,"Row and column permutations must be identity for in-place ILU"); 1637 } 1638 1639 outA = inA; 1640 inA->factor = FACTOR_LU; 1641 a->row = row; 1642 a->col = col; 1643 1644 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1645 ierr = ISInvertPermutation(col,&(a->icol));CHKERRQ(ierr); 1646 PLogObjectParent(inA,a->icol); 1647 1648 if (!a->solve_work) { /* this matrix may have been factored before */ 1649 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1650 } 1651 1652 if (!a->diag) { 1653 ierr = MatMarkDiag_SeqAIJ(inA);CHKERRQ(ierr); 1654 } 1655 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #include "pinclude/blaslapack.h" 1660 #undef __FUNC__ 1661 #define __FUNC__ "MatScale_SeqAIJ" 1662 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1663 { 1664 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1665 int one = 1; 1666 1667 PetscFunctionBegin; 1668 BLscal_( &a->nz, alpha, a->a, &one ); 1669 PLogFlops(a->nz); 1670 PetscFunctionReturn(0); 1671 } 1672 1673 #undef __FUNC__ 1674 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1675 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B) 1676 { 1677 int ierr,i; 1678 1679 PetscFunctionBegin; 1680 if (scall == MAT_INITIAL_MATRIX) { 1681 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) );CHKPTRQ(*B); 1682 } 1683 1684 for ( i=0; i<n; i++ ) { 1685 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1686 } 1687 PetscFunctionReturn(0); 1688 } 1689 1690 #undef __FUNC__ 1691 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1692 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1693 { 1694 PetscFunctionBegin; 1695 *bs = 1; 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNC__ 1700 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1701 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1702 { 1703 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1704 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1705 int start, end, *ai, *aj; 1706 BT table; 1707 1708 PetscFunctionBegin; 1709 shift = a->indexshift; 1710 m = a->m; 1711 ai = a->i; 1712 aj = a->j+shift; 1713 1714 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1715 1716 nidx = (int *) PetscMalloc((m+1)*sizeof(int));CHKPTRQ(nidx); 1717 ierr = BTCreate(m,table);CHKERRQ(ierr); 1718 1719 for ( i=0; i<is_max; i++ ) { 1720 /* Initialize the two local arrays */ 1721 isz = 0; 1722 ierr = BTMemzero(m,table);CHKERRQ(ierr); 1723 1724 /* Extract the indices, assume there can be duplicate entries */ 1725 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1726 ierr = ISGetSize(is[i],&n);CHKERRQ(ierr); 1727 1728 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1729 for ( j=0; j<n ; ++j){ 1730 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1731 } 1732 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1733 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1734 1735 k = 0; 1736 for ( j=0; j<ov; j++){ /* for each overlap */ 1737 n = isz; 1738 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1739 row = nidx[k]; 1740 start = ai[row]; 1741 end = ai[row+1]; 1742 for ( l = start; l<end ; l++){ 1743 val = aj[l] + shift; 1744 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1745 } 1746 } 1747 } 1748 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i));CHKERRQ(ierr); 1749 } 1750 ierr = BTDestroy(table);CHKERRQ(ierr); 1751 ierr = PetscFree(nidx);CHKERRQ(ierr); 1752 PetscFunctionReturn(0); 1753 } 1754 1755 /* -------------------------------------------------------------- */ 1756 #undef __FUNC__ 1757 #define __FUNC__ "MatPermute_SeqAIJ" 1758 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1759 { 1760 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1761 Scalar *vwork; 1762 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1763 int *row,*col,*cnew,j,*lens; 1764 IS icolp,irowp; 1765 1766 PetscFunctionBegin; 1767 ierr = ISInvertPermutation(rowp,&irowp);CHKERRQ(ierr); 1768 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1769 ierr = ISInvertPermutation(colp,&icolp);CHKERRQ(ierr); 1770 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1771 1772 /* determine lengths of permuted rows */ 1773 lens = (int *) PetscMalloc( (m+1)*sizeof(int) );CHKPTRQ(lens); 1774 for (i=0; i<m; i++ ) { 1775 lens[row[i]] = a->i[i+1] - a->i[i]; 1776 } 1777 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1778 ierr = PetscFree(lens);CHKERRQ(ierr); 1779 1780 cnew = (int *) PetscMalloc( n*sizeof(int) );CHKPTRQ(cnew); 1781 for (i=0; i<m; i++) { 1782 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1783 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1784 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1785 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1786 } 1787 ierr = PetscFree(cnew);CHKERRQ(ierr); 1788 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1789 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1790 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1791 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1792 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1793 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1794 PetscFunctionReturn(0); 1795 } 1796 1797 #undef __FUNC__ 1798 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1799 int MatPrintHelp_SeqAIJ(Mat A) 1800 { 1801 static int called = 0; 1802 MPI_Comm comm = A->comm; 1803 int ierr; 1804 1805 PetscFunctionBegin; 1806 if (called) {PetscFunctionReturn(0);} else called = 1; 1807 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1808 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1809 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1810 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1811 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1812 #if defined(PETSC_HAVE_ESSL) 1813 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1814 #endif 1815 PetscFunctionReturn(0); 1816 } 1817 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1818 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1819 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1820 1821 #undef __FUNC__ 1822 #define __FUNC__ "MatCopy_SeqAIJ" 1823 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1824 { 1825 int ierr; 1826 1827 PetscFunctionBegin; 1828 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1829 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1830 Mat_SeqAIJ *b = (Mat_SeqAIJ *) B->data; 1831 1832 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1833 SETERRQ(1,1,"Number of nonzeros in two matrices are different"); 1834 } 1835 ierr = PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 1836 } else { 1837 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1838 } 1839 PetscFunctionReturn(0); 1840 } 1841 1842 /* -------------------------------------------------------------------*/ 1843 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1844 MatGetRow_SeqAIJ, 1845 MatRestoreRow_SeqAIJ, 1846 MatMult_SeqAIJ, 1847 MatMultAdd_SeqAIJ, 1848 MatMultTrans_SeqAIJ, 1849 MatMultTransAdd_SeqAIJ, 1850 MatSolve_SeqAIJ, 1851 MatSolveAdd_SeqAIJ, 1852 MatSolveTrans_SeqAIJ, 1853 MatSolveTransAdd_SeqAIJ, 1854 MatLUFactor_SeqAIJ, 1855 0, 1856 MatRelax_SeqAIJ, 1857 MatTranspose_SeqAIJ, 1858 MatGetInfo_SeqAIJ, 1859 MatEqual_SeqAIJ, 1860 MatGetDiagonal_SeqAIJ, 1861 MatDiagonalScale_SeqAIJ, 1862 MatNorm_SeqAIJ, 1863 0, 1864 MatAssemblyEnd_SeqAIJ, 1865 MatCompress_SeqAIJ, 1866 MatSetOption_SeqAIJ, 1867 MatZeroEntries_SeqAIJ, 1868 MatZeroRows_SeqAIJ, 1869 MatLUFactorSymbolic_SeqAIJ, 1870 MatLUFactorNumeric_SeqAIJ, 1871 0, 1872 0, 1873 MatGetSize_SeqAIJ, 1874 MatGetSize_SeqAIJ, 1875 MatGetOwnershipRange_SeqAIJ, 1876 MatILUFactorSymbolic_SeqAIJ, 1877 0, 1878 0, 1879 0, 1880 MatDuplicate_SeqAIJ, 1881 0, 1882 0, 1883 MatILUFactor_SeqAIJ, 1884 0, 1885 0, 1886 MatGetSubMatrices_SeqAIJ, 1887 MatIncreaseOverlap_SeqAIJ, 1888 MatGetValues_SeqAIJ, 1889 MatCopy_SeqAIJ, 1890 MatPrintHelp_SeqAIJ, 1891 MatScale_SeqAIJ, 1892 0, 1893 0, 1894 MatILUDTFactor_SeqAIJ, 1895 MatGetBlockSize_SeqAIJ, 1896 MatGetRowIJ_SeqAIJ, 1897 MatRestoreRowIJ_SeqAIJ, 1898 MatGetColumnIJ_SeqAIJ, 1899 MatRestoreColumnIJ_SeqAIJ, 1900 MatFDColoringCreate_SeqAIJ, 1901 MatColoringPatch_SeqAIJ, 1902 0, 1903 MatPermute_SeqAIJ, 1904 0, 1905 0, 1906 0, 1907 0, 1908 MatGetMaps_Petsc}; 1909 1910 extern int MatUseSuperLU_SeqAIJ(Mat); 1911 extern int MatUseEssl_SeqAIJ(Mat); 1912 extern int MatUseDXML_SeqAIJ(Mat); 1913 1914 EXTERN_C_BEGIN 1915 #undef __FUNC__ 1916 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1917 1918 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1919 { 1920 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1921 int i,nz,n; 1922 1923 PetscFunctionBegin; 1924 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1925 1926 nz = aij->maxnz; 1927 n = aij->n; 1928 for (i=0; i<nz; i++) { 1929 aij->j[i] = indices[i]; 1930 } 1931 aij->nz = nz; 1932 for ( i=0; i<n; i++ ) { 1933 aij->ilen[i] = aij->imax[i]; 1934 } 1935 1936 PetscFunctionReturn(0); 1937 } 1938 EXTERN_C_END 1939 1940 #undef __FUNC__ 1941 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1942 /*@ 1943 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1944 in the matrix. 1945 1946 Input Parameters: 1947 + mat - the SeqAIJ matrix 1948 - indices - the column indices 1949 1950 Level: advanced 1951 1952 Notes: 1953 This can be called if you have precomputed the nonzero structure of the 1954 matrix and want to provide it to the matrix object to improve the performance 1955 of the MatSetValues() operation. 1956 1957 You MUST have set the correct numbers of nonzeros per row in the call to 1958 MatCreateSeqAIJ(). 1959 1960 MUST be called before any calls to MatSetValues(); 1961 1962 @*/ 1963 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1964 { 1965 int ierr,(*f)(Mat,int *); 1966 1967 PetscFunctionBegin; 1968 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1969 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1970 if (f) { 1971 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1972 } else { 1973 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1974 } 1975 PetscFunctionReturn(0); 1976 } 1977 1978 /* ----------------------------------------------------------------------------------------*/ 1979 1980 EXTERN_C_BEGIN 1981 #undef __FUNC__ 1982 #define __FUNC__ "MatStoreValues_SeqAIJ" 1983 int MatStoreValues_SeqAIJ(Mat mat) 1984 { 1985 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1986 int nz = aij->i[aij->m]+aij->indexshift,ierr; 1987 1988 PetscFunctionBegin; 1989 if (aij->nonew != 1) { 1990 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1991 } 1992 1993 /* allocate space for values if not already there */ 1994 if (!aij->saved_values) { 1995 aij->saved_values = (Scalar *) PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1996 } 1997 1998 /* copy values over */ 1999 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 2000 PetscFunctionReturn(0); 2001 } 2002 EXTERN_C_END 2003 2004 #undef __FUNC__ 2005 #define __FUNC__ "MatStoreValues" 2006 /*@ 2007 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2008 example, reuse of the linear part of a Jacobian, while recomputing the 2009 nonlinear portion. 2010 2011 Collect on Mat 2012 2013 Input Parameters: 2014 . mat - the matrix (currently on AIJ matrices support this option) 2015 2016 Level: advanced 2017 2018 Common Usage, with SNESSolve(): 2019 $ Create Jacobian matrix 2020 $ Set linear terms into matrix 2021 $ Apply boundary conditions to matrix, at this time matrix must have 2022 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2023 $ boundary conditions again will not change the nonzero structure 2024 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2025 $ ierr = MatStoreValues(mat); 2026 $ Call SNESSetJacobian() with matrix 2027 $ In your Jacobian routine 2028 $ ierr = MatRetrieveValues(mat); 2029 $ Set nonlinear terms in matrix 2030 2031 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2032 $ // build linear portion of Jacobian 2033 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2034 $ ierr = MatStoreValues(mat); 2035 $ loop over nonlinear iterations 2036 $ ierr = MatRetrieveValues(mat); 2037 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2038 $ // call MatAssemblyBegin/End() on matrix 2039 $ Solve linear system with Jacobian 2040 $ endloop 2041 2042 Notes: 2043 Matrix must already be assemblied before calling this routine 2044 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2045 calling this routine. 2046 2047 .seealso: MatRetrieveValues() 2048 2049 @*/ 2050 int MatStoreValues(Mat mat) 2051 { 2052 int ierr,(*f)(Mat); 2053 2054 PetscFunctionBegin; 2055 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2056 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2057 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2058 2059 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f);CHKERRQ(ierr); 2060 if (f) { 2061 ierr = (*f)(mat);CHKERRQ(ierr); 2062 } else { 2063 SETERRQ(1,1,"Wrong type of matrix to store values"); 2064 } 2065 PetscFunctionReturn(0); 2066 } 2067 2068 EXTERN_C_BEGIN 2069 #undef __FUNC__ 2070 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 2071 int MatRetrieveValues_SeqAIJ(Mat mat) 2072 { 2073 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2074 int nz = aij->i[aij->m]+aij->indexshift,ierr; 2075 2076 PetscFunctionBegin; 2077 if (aij->nonew != 1) { 2078 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2079 } 2080 if (!aij->saved_values) { 2081 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 2082 } 2083 2084 /* copy values over */ 2085 ierr = PetscMemcpy(aij->a, aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 2086 PetscFunctionReturn(0); 2087 } 2088 EXTERN_C_END 2089 2090 #undef __FUNC__ 2091 #define __FUNC__ "MatRetrieveValues" 2092 /*@ 2093 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2094 example, reuse of the linear part of a Jacobian, while recomputing the 2095 nonlinear portion. 2096 2097 Collect on Mat 2098 2099 Input Parameters: 2100 . mat - the matrix (currently on AIJ matrices support this option) 2101 2102 Level: advanced 2103 2104 .seealso: MatStoreValues() 2105 2106 @*/ 2107 int MatRetrieveValues(Mat mat) 2108 { 2109 int ierr,(*f)(Mat); 2110 2111 PetscFunctionBegin; 2112 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2113 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2114 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2115 2116 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f);CHKERRQ(ierr); 2117 if (f) { 2118 ierr = (*f)(mat);CHKERRQ(ierr); 2119 } else { 2120 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 2121 } 2122 PetscFunctionReturn(0); 2123 } 2124 2125 /* --------------------------------------------------------------------------------*/ 2126 2127 #undef __FUNC__ 2128 #define __FUNC__ "MatCreateSeqAIJ" 2129 /*@C 2130 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2131 (the default parallel PETSc format). For good matrix assembly performance 2132 the user should preallocate the matrix storage by setting the parameter nz 2133 (or the array nnz). By setting these parameters accurately, performance 2134 during matrix assembly can be increased by more than a factor of 50. 2135 2136 Collective on MPI_Comm 2137 2138 Input Parameters: 2139 + comm - MPI communicator, set to PETSC_COMM_SELF 2140 . m - number of rows 2141 . n - number of columns 2142 . nz - number of nonzeros per row (same for all rows) 2143 - nnz - array containing the number of nonzeros in the various rows 2144 (possibly different for each row) or PETSC_NULL 2145 2146 Output Parameter: 2147 . A - the matrix 2148 2149 Notes: 2150 The AIJ format (also called the Yale sparse matrix format or 2151 compressed row storage), is fully compatible with standard Fortran 77 2152 storage. That is, the stored row and column indices can begin at 2153 either one (as in Fortran) or zero. See the users' manual for details. 2154 2155 Specify the preallocated storage with either nz or nnz (not both). 2156 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2157 allocation. For large problems you MUST preallocate memory or you 2158 will get TERRIBLE performance, see the users' manual chapter on matrices. 2159 2160 By default, this format uses inodes (identical nodes) when possible, to 2161 improve numerical efficiency of matrix-vector products and solves. We 2162 search for consecutive rows with the same nonzero structure, thereby 2163 reusing matrix information to achieve increased efficiency. 2164 2165 Options Database Keys: 2166 + -mat_aij_no_inode - Do not use inodes 2167 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2168 - -mat_aij_oneindex - Internally use indexing starting at 1 2169 rather than 0. Note that when calling MatSetValues(), 2170 the user still MUST index entries starting at 0! 2171 2172 Level: intermediate 2173 2174 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 2175 @*/ 2176 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 2177 { 2178 Mat B; 2179 Mat_SeqAIJ *b; 2180 int i, len, ierr, flg,size; 2181 2182 PetscFunctionBegin; 2183 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2184 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2185 2186 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"nz cannot be less than -2: value %d",nz); 2187 if (nnz) { 2188 for (i=0; i<m; i++) { 2189 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,0,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2190 } 2191 } 2192 2193 *A = 0; 2194 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",comm,MatDestroy,MatView); 2195 PLogObjectCreate(B); 2196 B->data = (void *) (b = PetscNew(Mat_SeqAIJ));CHKPTRQ(b); 2197 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2198 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2199 B->ops->destroy = MatDestroy_SeqAIJ; 2200 B->ops->view = MatView_SeqAIJ; 2201 B->factor = 0; 2202 B->lupivotthreshold = 1.0; 2203 B->mapping = 0; 2204 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 2205 b->ilu_preserve_row_sums = PETSC_FALSE; 2206 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2207 b->row = 0; 2208 b->col = 0; 2209 b->icol = 0; 2210 b->indexshift = 0; 2211 b->reallocs = 0; 2212 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg);CHKERRQ(ierr); 2213 if (flg) b->indexshift = -1; 2214 2215 b->m = m; B->m = m; B->M = m; 2216 b->n = n; B->n = n; B->N = n; 2217 2218 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2219 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2220 2221 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) );CHKPTRQ(b->imax); 2222 if (nnz == PETSC_NULL) { 2223 if (nz == PETSC_DEFAULT) nz = 10; 2224 else if (nz <= 0) nz = 1; 2225 for ( i=0; i<m; i++ ) b->imax[i] = nz; 2226 nz = nz*m; 2227 } else { 2228 nz = 0; 2229 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 2230 } 2231 2232 /* allocate the matrix space */ 2233 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2234 b->a = (Scalar *) PetscMalloc( len );CHKPTRQ(b->a); 2235 b->j = (int *) (b->a + nz); 2236 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2237 b->i = b->j + nz; 2238 b->singlemalloc = 1; 2239 2240 b->i[0] = -b->indexshift; 2241 for (i=1; i<m+1; i++) { 2242 b->i[i] = b->i[i-1] + b->imax[i-1]; 2243 } 2244 2245 /* b->ilen will count nonzeros in each row so far. */ 2246 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 2247 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2248 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 2249 2250 b->nz = 0; 2251 b->maxnz = nz; 2252 b->sorted = 0; 2253 b->roworiented = 1; 2254 b->nonew = 0; 2255 b->diag = 0; 2256 b->solve_work = 0; 2257 b->spptr = 0; 2258 b->inode.node_count = 0; 2259 b->inode.size = 0; 2260 b->inode.limit = 5; 2261 b->inode.max_limit = 5; 2262 b->saved_values = 0; 2263 B->info.nz_unneeded = (double)b->maxnz; 2264 b->idiag = 0; 2265 b->ssor = 0; 2266 2267 *A = B; 2268 2269 /* SuperLU is not currently supported through PETSc */ 2270 #if defined(PETSC_HAVE_SUPERLU) 2271 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg);CHKERRQ(ierr); 2272 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2273 #endif 2274 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg);CHKERRQ(ierr); 2275 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2276 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg);CHKERRQ(ierr); 2277 if (flg) { 2278 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2279 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2280 } 2281 ierr = OptionsHasName(PETSC_NULL,"-help", &flg);CHKERRQ(ierr); 2282 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr); } 2283 2284 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2285 "MatSeqAIJSetColumnIndices_SeqAIJ", 2286 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C", 2288 "MatStoreValues_SeqAIJ", 2289 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C", 2291 "MatRetrieveValues_SeqAIJ", 2292 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2293 PetscFunctionReturn(0); 2294 } 2295 2296 #undef __FUNC__ 2297 #define __FUNC__ "MatDuplicate_SeqAIJ" 2298 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2299 { 2300 Mat C; 2301 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 2302 int i,len, m = a->m,shift = a->indexshift,ierr; 2303 2304 PetscFunctionBegin; 2305 *B = 0; 2306 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",A->comm,MatDestroy,MatView); 2307 PLogObjectCreate(C); 2308 C->data = (void *) (c = PetscNew(Mat_SeqAIJ));CHKPTRQ(c); 2309 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2310 C->ops->destroy = MatDestroy_SeqAIJ; 2311 C->ops->view = MatView_SeqAIJ; 2312 C->factor = A->factor; 2313 c->row = 0; 2314 c->col = 0; 2315 c->icol = 0; 2316 c->indexshift = shift; 2317 C->assembled = PETSC_TRUE; 2318 2319 c->m = C->m = a->m; 2320 c->n = C->n = a->n; 2321 C->M = a->m; 2322 C->N = a->n; 2323 2324 c->imax = (int *) PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->imax); 2325 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->ilen); 2326 for ( i=0; i<m; i++ ) { 2327 c->imax[i] = a->imax[i]; 2328 c->ilen[i] = a->ilen[i]; 2329 } 2330 2331 /* allocate the matrix space */ 2332 c->singlemalloc = 1; 2333 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2334 c->a = (Scalar *) PetscMalloc( len );CHKPTRQ(c->a); 2335 c->j = (int *) (c->a + a->i[m] + shift); 2336 c->i = c->j + a->i[m] + shift; 2337 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2338 if (m > 0) { 2339 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2340 if (cpvalues == MAT_COPY_VALUES) { 2341 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2342 } else { 2343 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2344 } 2345 } 2346 2347 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2348 c->sorted = a->sorted; 2349 c->roworiented = a->roworiented; 2350 c->nonew = a->nonew; 2351 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2352 c->saved_values = 0; 2353 c->idiag = 0; 2354 c->ssor = 0; 2355 2356 if (a->diag) { 2357 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) );CHKPTRQ(c->diag); 2358 PLogObjectMemory(C,(m+1)*sizeof(int)); 2359 for ( i=0; i<m; i++ ) { 2360 c->diag[i] = a->diag[i]; 2361 } 2362 } else c->diag = 0; 2363 c->inode.limit = a->inode.limit; 2364 c->inode.max_limit = a->inode.max_limit; 2365 if (a->inode.size){ 2366 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) );CHKPTRQ(c->inode.size); 2367 c->inode.node_count = a->inode.node_count; 2368 ierr = PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int));CHKERRQ(ierr); 2369 } else { 2370 c->inode.size = 0; 2371 c->inode.node_count = 0; 2372 } 2373 c->nz = a->nz; 2374 c->maxnz = a->maxnz; 2375 c->solve_work = 0; 2376 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2377 2378 *B = C; 2379 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2380 PetscFunctionReturn(0); 2381 } 2382 2383 #undef __FUNC__ 2384 #define __FUNC__ "MatLoad_SeqAIJ" 2385 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2386 { 2387 Mat_SeqAIJ *a; 2388 Mat B; 2389 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2390 MPI_Comm comm; 2391 2392 PetscFunctionBegin; 2393 ierr = PetscObjectGetComm((PetscObject) viewer,&comm);CHKERRQ(ierr); 2394 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2395 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2396 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2397 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2398 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2399 M = header[1]; N = header[2]; nz = header[3]; 2400 2401 if (nz < 0) { 2402 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2403 } 2404 2405 /* read in row lengths */ 2406 rowlengths = (int*) PetscMalloc( M*sizeof(int) );CHKPTRQ(rowlengths); 2407 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2408 2409 /* create our matrix */ 2410 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2411 B = *A; 2412 a = (Mat_SeqAIJ *) B->data; 2413 shift = a->indexshift; 2414 2415 /* read in column indices and adjust for Fortran indexing*/ 2416 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2417 if (shift) { 2418 for ( i=0; i<nz; i++ ) { 2419 a->j[i] += 1; 2420 } 2421 } 2422 2423 /* read in nonzero values */ 2424 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2425 2426 /* set matrix "i" values */ 2427 a->i[0] = -shift; 2428 for ( i=1; i<= M; i++ ) { 2429 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2430 a->ilen[i-1] = rowlengths[i-1]; 2431 } 2432 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2433 2434 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2435 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2436 PetscFunctionReturn(0); 2437 } 2438 2439 #undef __FUNC__ 2440 #define __FUNC__ "MatEqual_SeqAIJ" 2441 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2442 { 2443 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2444 2445 PetscFunctionBegin; 2446 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2447 2448 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2449 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2450 (a->indexshift != b->indexshift)) { 2451 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2452 } 2453 2454 /* if the a->i are the same */ 2455 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2456 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2457 } 2458 2459 /* if a->j are the same */ 2460 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2461 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2462 } 2463 2464 /* if a->a are the same */ 2465 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2466 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2467 } 2468 *flg = PETSC_TRUE; 2469 PetscFunctionReturn(0); 2470 2471 } 2472