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