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