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