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