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