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