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