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