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