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