1 /*$Id: aij.c,v 1.355 2000/09/22 20:43:38 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,"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,"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,"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,"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,"Negative row: %d",row); 239 if (row >= a->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"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,"Negative column: %d",in[l]); 244 if (in[l] >= a->n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"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,"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,"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,"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,"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,"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,"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,"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,MatLUInfo*,Mat*); 1223 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1224 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 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 = ISGetLocalSize(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,"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,"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,"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,"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,"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,"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; 1409 aj += len; 1410 } 1411 if (shift) { 1412 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1413 } 1414 1415 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1416 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1417 1418 if (B) { 1419 *B = C; 1420 } else { 1421 PetscOps *Abops; 1422 MatOps Aops; 1423 1424 /* This isn't really an in-place transpose */ 1425 ierr = PetscFree(a->a);CHKERRQ(ierr); 1426 if (!a->singlemalloc) { 1427 ierr = PetscFree(a->i);CHKERRQ(ierr); 1428 ierr = PetscFree(a->j); 1429 } 1430 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 1431 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 1432 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 1433 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 1434 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 1435 ierr = PetscFree(a);CHKERRQ(ierr); 1436 1437 1438 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1439 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1440 1441 /* 1442 This is horrible, horrible code. We need to keep the 1443 the bops and ops Structures, copy everything from C 1444 including the function pointers.. 1445 */ 1446 ierr = PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1447 ierr = PetscMemcpy(A->bops,C->bops,sizeof(PetscOps));CHKERRQ(ierr); 1448 Abops = A->bops; 1449 Aops = A->ops; 1450 refct = A->refct; 1451 ierr = PetscMemcpy(A,C,sizeof(struct _p_Mat));CHKERRQ(ierr); 1452 A->bops = Abops; 1453 A->ops = Aops; 1454 A->qlist = 0; 1455 A->refct = refct; 1456 /* copy over the type_name and name */ 1457 ierr = PetscStrallocpy(C->type_name,&A->type_name);CHKERRQ(ierr); 1458 ierr = PetscStrallocpy(C->name,&A->name);CHKERRQ(ierr); 1459 1460 PetscHeaderDestroy(C); 1461 } 1462 PetscFunctionReturn(0); 1463 } 1464 1465 #undef __FUNC__ 1466 #define __FUNC__ /*<a name="MatDiagonalScale_SeqAIJ"></a>*/"MatDiagonalScale_SeqAIJ" 1467 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1468 { 1469 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1470 Scalar *l,*r,x,*v; 1471 int ierr,i,j,m = a->m,n = a->n,M,nz = a->nz,*jj,shift = a->indexshift; 1472 1473 PetscFunctionBegin; 1474 if (ll) { 1475 /* The local size is used so that VecMPI can be passed to this routine 1476 by MatDiagonalScale_MPIAIJ */ 1477 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1478 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1479 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1480 v = a->a; 1481 for (i=0; i<m; i++) { 1482 x = l[i]; 1483 M = a->i[i+1] - a->i[i]; 1484 for (j=0; j<M; j++) { (*v++) *= x;} 1485 } 1486 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1487 PLogFlops(nz); 1488 } 1489 if (rr) { 1490 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1491 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1492 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1493 v = a->a; jj = a->j; 1494 for (i=0; i<nz; i++) { 1495 (*v++) *= r[*jj++ + shift]; 1496 } 1497 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1498 PLogFlops(nz); 1499 } 1500 PetscFunctionReturn(0); 1501 } 1502 1503 #undef __FUNC__ 1504 #define __FUNC__ /*<a name="MatGetSubMatrix_SeqAIJ"></a>*/"MatGetSubMatrix_SeqAIJ" 1505 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1506 { 1507 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1508 int *smap,i,k,kstart,kend,ierr,oldcols = a->n,*lens; 1509 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1510 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1511 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1512 Scalar *a_new,*mat_a; 1513 Mat C; 1514 PetscTruth stride; 1515 1516 PetscFunctionBegin; 1517 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1518 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1519 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1520 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1521 1522 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1523 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1524 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1525 1526 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1527 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1528 if (stride && step == 1) { 1529 /* special case of contiguous rows */ 1530 lens = (int*)PetscMalloc((ncols+nrows+1)*sizeof(int));CHKPTRQ(lens); 1531 starts = lens + ncols; 1532 /* loop over new rows determining lens and starting points */ 1533 for (i=0; i<nrows; i++) { 1534 kstart = ai[irow[i]]+shift; 1535 kend = kstart + ailen[irow[i]]; 1536 for (k=kstart; k<kend; k++) { 1537 if (aj[k]+shift >= first) { 1538 starts[i] = k; 1539 break; 1540 } 1541 } 1542 sum = 0; 1543 while (k < kend) { 1544 if (aj[k++]+shift >= first+ncols) break; 1545 sum++; 1546 } 1547 lens[i] = sum; 1548 } 1549 /* create submatrix */ 1550 if (scall == MAT_REUSE_MATRIX) { 1551 int n_cols,n_rows; 1552 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1553 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1554 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1555 C = *B; 1556 } else { 1557 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1558 } 1559 c = (Mat_SeqAIJ*)C->data; 1560 1561 /* loop over rows inserting into submatrix */ 1562 a_new = c->a; 1563 j_new = c->j; 1564 i_new = c->i; 1565 i_new[0] = -shift; 1566 for (i=0; i<nrows; i++) { 1567 ii = starts[i]; 1568 lensi = lens[i]; 1569 for (k=0; k<lensi; k++) { 1570 *j_new++ = aj[ii+k] - first; 1571 } 1572 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar));CHKERRQ(ierr); 1573 a_new += lensi; 1574 i_new[i+1] = i_new[i] + lensi; 1575 c->ilen[i] = lensi; 1576 } 1577 ierr = PetscFree(lens);CHKERRQ(ierr); 1578 } else { 1579 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1580 smap = (int*)PetscMalloc((1+oldcols)*sizeof(int));CHKPTRQ(smap); 1581 ssmap = smap + shift; 1582 lens = (int*)PetscMalloc((1+nrows)*sizeof(int));CHKPTRQ(lens); 1583 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1584 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1585 /* determine lens of each row */ 1586 for (i=0; i<nrows; i++) { 1587 kstart = ai[irow[i]]+shift; 1588 kend = kstart + a->ilen[irow[i]]; 1589 lens[i] = 0; 1590 for (k=kstart; k<kend; k++) { 1591 if (ssmap[aj[k]]) { 1592 lens[i]++; 1593 } 1594 } 1595 } 1596 /* Create and fill new matrix */ 1597 if (scall == MAT_REUSE_MATRIX) { 1598 PetscTruth equal; 1599 1600 c = (Mat_SeqAIJ *)((*B)->data); 1601 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1602 ierr = PetscMemcmp(c->ilen,lens,c->m*sizeof(int),&equal);CHKERRQ(ierr); 1603 if (!equal) { 1604 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1605 } 1606 ierr = PetscMemzero(c->ilen,c->m*sizeof(int));CHKERRQ(ierr); 1607 C = *B; 1608 } else { 1609 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1610 } 1611 c = (Mat_SeqAIJ *)(C->data); 1612 for (i=0; i<nrows; i++) { 1613 row = irow[i]; 1614 kstart = ai[row]+shift; 1615 kend = kstart + a->ilen[row]; 1616 mat_i = c->i[i]+shift; 1617 mat_j = c->j + mat_i; 1618 mat_a = c->a + mat_i; 1619 mat_ilen = c->ilen + i; 1620 for (k=kstart; k<kend; k++) { 1621 if ((tcol=ssmap[a->j[k]])) { 1622 *mat_j++ = tcol - (!shift); 1623 *mat_a++ = a->a[k]; 1624 (*mat_ilen)++; 1625 1626 } 1627 } 1628 } 1629 /* Free work space */ 1630 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1631 ierr = PetscFree(smap);CHKERRQ(ierr); 1632 ierr = PetscFree(lens);CHKERRQ(ierr); 1633 } 1634 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1635 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1636 1637 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1638 *B = C; 1639 PetscFunctionReturn(0); 1640 } 1641 1642 /* 1643 */ 1644 #undef __FUNC__ 1645 #define __FUNC__ /*<a name="MatILUFactor_SeqAIJ"></a>*/"MatILUFactor_SeqAIJ" 1646 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1647 { 1648 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1649 int ierr; 1650 Mat outA; 1651 PetscTruth row_identity,col_identity; 1652 1653 PetscFunctionBegin; 1654 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1655 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1656 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1657 if (!row_identity || !col_identity) { 1658 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1659 } 1660 1661 outA = inA; 1662 inA->factor = FACTOR_LU; 1663 a->row = row; 1664 a->col = col; 1665 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1666 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1667 1668 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1669 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1670 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1671 PLogObjectParent(inA,a->icol); 1672 1673 if (!a->solve_work) { /* this matrix may have been factored before */ 1674 a->solve_work = (Scalar*)PetscMalloc((a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1675 } 1676 1677 if (!a->diag) { 1678 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1679 } 1680 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #include "petscblaslapack.h" 1685 #undef __FUNC__ 1686 #define __FUNC__ /*<a name="MatScale_SeqAIJ"></a>*/"MatScale_SeqAIJ" 1687 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1688 { 1689 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1690 int one = 1; 1691 1692 PetscFunctionBegin; 1693 BLscal_(&a->nz,alpha,a->a,&one); 1694 PLogFlops(a->nz); 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNC__ 1699 #define __FUNC__ /*<a name="MatGetSubMatrices_SeqAIJ"></a>*/"MatGetSubMatrices_SeqAIJ" 1700 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1701 { 1702 int ierr,i; 1703 1704 PetscFunctionBegin; 1705 if (scall == MAT_INITIAL_MATRIX) { 1706 *B = (Mat*)PetscMalloc((n+1)*sizeof(Mat));CHKPTRQ(*B); 1707 } 1708 1709 for (i=0; i<n; i++) { 1710 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1711 } 1712 PetscFunctionReturn(0); 1713 } 1714 1715 #undef __FUNC__ 1716 #define __FUNC__ /*<a name="MatGetBlockSize_SeqAIJ"></a>*/"MatGetBlockSize_SeqAIJ" 1717 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1718 { 1719 PetscFunctionBegin; 1720 *bs = 1; 1721 PetscFunctionReturn(0); 1722 } 1723 1724 #undef __FUNC__ 1725 #define __FUNC__ /*<a name="MatIncreaseOverlap_SeqAIJ"></a>*/"MatIncreaseOverlap_SeqAIJ" 1726 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1727 { 1728 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1729 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1730 int start,end,*ai,*aj; 1731 PetscBT table; 1732 1733 PetscFunctionBegin; 1734 shift = a->indexshift; 1735 m = a->m; 1736 ai = a->i; 1737 aj = a->j+shift; 1738 1739 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1740 1741 nidx = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(nidx); 1742 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1743 1744 for (i=0; i<is_max; i++) { 1745 /* Initialize the two local arrays */ 1746 isz = 0; 1747 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1748 1749 /* Extract the indices, assume there can be duplicate entries */ 1750 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1751 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1752 1753 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1754 for (j=0; j<n ; ++j){ 1755 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1756 } 1757 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1758 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1759 1760 k = 0; 1761 for (j=0; j<ov; j++){ /* for each overlap */ 1762 n = isz; 1763 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1764 row = nidx[k]; 1765 start = ai[row]; 1766 end = ai[row+1]; 1767 for (l = start; l<end ; l++){ 1768 val = aj[l] + shift; 1769 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1770 } 1771 } 1772 } 1773 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1774 } 1775 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1776 ierr = PetscFree(nidx);CHKERRQ(ierr); 1777 PetscFunctionReturn(0); 1778 } 1779 1780 /* -------------------------------------------------------------- */ 1781 #undef __FUNC__ 1782 #define __FUNC__ /*<a name="MatPermute_SeqAIJ"></a>*/"MatPermute_SeqAIJ" 1783 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1784 { 1785 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1786 Scalar *vwork; 1787 int i,ierr,nz,m = a->m,n = a->n,*cwork; 1788 int *row,*col,*cnew,j,*lens; 1789 IS icolp,irowp; 1790 1791 PetscFunctionBegin; 1792 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1793 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1794 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1795 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1796 1797 /* determine lengths of permuted rows */ 1798 lens = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(lens); 1799 for (i=0; i<m; i++) { 1800 lens[row[i]] = a->i[i+1] - a->i[i]; 1801 } 1802 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1803 ierr = PetscFree(lens);CHKERRQ(ierr); 1804 1805 cnew = (int*)PetscMalloc(n*sizeof(int));CHKPTRQ(cnew); 1806 for (i=0; i<m; i++) { 1807 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1808 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1809 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1810 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1811 } 1812 ierr = PetscFree(cnew);CHKERRQ(ierr); 1813 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1814 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1815 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1816 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1817 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1818 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1819 PetscFunctionReturn(0); 1820 } 1821 1822 #undef __FUNC__ 1823 #define __FUNC__ /*<a name="MatPrintHelp_SeqAIJ"></a>*/"MatPrintHelp_SeqAIJ" 1824 int MatPrintHelp_SeqAIJ(Mat A) 1825 { 1826 static PetscTruth called = PETSC_FALSE; 1827 MPI_Comm comm = A->comm; 1828 int ierr; 1829 1830 PetscFunctionBegin; 1831 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1832 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1833 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1834 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1835 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1836 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1837 #if defined(PETSC_HAVE_ESSL) 1838 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1839 #endif 1840 #if defined(PETSC_HAVE_LUSOL) 1841 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1842 #endif 1843 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) 1844 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1845 #endif 1846 PetscFunctionReturn(0); 1847 } 1848 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1849 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1850 EXTERN int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1851 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1852 #undef __FUNC__ 1853 #define __FUNC__ /*<a name="MatCopy_SeqAIJ"></a>*/"MatCopy_SeqAIJ" 1854 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1855 { 1856 int ierr; 1857 1858 PetscFunctionBegin; 1859 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1860 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1861 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1862 1863 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1864 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1865 } 1866 ierr = PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 1867 } else { 1868 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /* -------------------------------------------------------------------*/ 1874 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1875 MatGetRow_SeqAIJ, 1876 MatRestoreRow_SeqAIJ, 1877 MatMult_SeqAIJ, 1878 MatMultAdd_SeqAIJ, 1879 MatMultTranspose_SeqAIJ, 1880 MatMultTransposeAdd_SeqAIJ, 1881 MatSolve_SeqAIJ, 1882 MatSolveAdd_SeqAIJ, 1883 MatSolveTranspose_SeqAIJ, 1884 MatSolveTransposeAdd_SeqAIJ, 1885 MatLUFactor_SeqAIJ, 1886 0, 1887 MatRelax_SeqAIJ, 1888 MatTranspose_SeqAIJ, 1889 MatGetInfo_SeqAIJ, 1890 MatEqual_SeqAIJ, 1891 MatGetDiagonal_SeqAIJ, 1892 MatDiagonalScale_SeqAIJ, 1893 MatNorm_SeqAIJ, 1894 0, 1895 MatAssemblyEnd_SeqAIJ, 1896 MatCompress_SeqAIJ, 1897 MatSetOption_SeqAIJ, 1898 MatZeroEntries_SeqAIJ, 1899 MatZeroRows_SeqAIJ, 1900 MatLUFactorSymbolic_SeqAIJ, 1901 MatLUFactorNumeric_SeqAIJ, 1902 0, 1903 0, 1904 MatGetSize_SeqAIJ, 1905 MatGetSize_SeqAIJ, 1906 MatGetOwnershipRange_SeqAIJ, 1907 MatILUFactorSymbolic_SeqAIJ, 1908 0, 1909 0, 1910 0, 1911 MatDuplicate_SeqAIJ, 1912 0, 1913 0, 1914 MatILUFactor_SeqAIJ, 1915 0, 1916 0, 1917 MatGetSubMatrices_SeqAIJ, 1918 MatIncreaseOverlap_SeqAIJ, 1919 MatGetValues_SeqAIJ, 1920 MatCopy_SeqAIJ, 1921 MatPrintHelp_SeqAIJ, 1922 MatScale_SeqAIJ, 1923 0, 1924 0, 1925 MatILUDTFactor_SeqAIJ, 1926 MatGetBlockSize_SeqAIJ, 1927 MatGetRowIJ_SeqAIJ, 1928 MatRestoreRowIJ_SeqAIJ, 1929 MatGetColumnIJ_SeqAIJ, 1930 MatRestoreColumnIJ_SeqAIJ, 1931 MatFDColoringCreate_SeqAIJ, 1932 MatColoringPatch_SeqAIJ, 1933 0, 1934 MatPermute_SeqAIJ, 1935 0, 1936 0, 1937 MatDestroy_SeqAIJ, 1938 MatView_SeqAIJ, 1939 MatGetMaps_Petsc}; 1940 1941 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 1942 EXTERN int MatUseEssl_SeqAIJ(Mat); 1943 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 1944 EXTERN int MatUseMatlab_SeqAIJ(Mat); 1945 EXTERN int MatUseDXML_SeqAIJ(Mat); 1946 1947 EXTERN_C_BEGIN 1948 #undef __FUNC__ 1949 #define __FUNC__ /*<a name="MatSeqAIJSetColumnIndices_SeqAIJ"></a>*/"MatSeqAIJSetColumnIndices_SeqAIJ" 1950 1951 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1952 { 1953 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1954 int i,nz,n; 1955 1956 PetscFunctionBegin; 1957 if (aij->indexshift) SETERRQ(1,"No support with 1 based indexing"); 1958 1959 nz = aij->maxnz; 1960 n = aij->n; 1961 for (i=0; i<nz; i++) { 1962 aij->j[i] = indices[i]; 1963 } 1964 aij->nz = nz; 1965 for (i=0; i<n; i++) { 1966 aij->ilen[i] = aij->imax[i]; 1967 } 1968 1969 PetscFunctionReturn(0); 1970 } 1971 EXTERN_C_END 1972 1973 #undef __FUNC__ 1974 #define __FUNC__ /*<a name="MatSeqAIJSetColumnIndices"></a>*/"MatSeqAIJSetColumnIndices" 1975 /*@ 1976 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1977 in the matrix. 1978 1979 Input Parameters: 1980 + mat - the SeqAIJ matrix 1981 - indices - the column indices 1982 1983 Level: advanced 1984 1985 Notes: 1986 This can be called if you have precomputed the nonzero structure of the 1987 matrix and want to provide it to the matrix object to improve the performance 1988 of the MatSetValues() operation. 1989 1990 You MUST have set the correct numbers of nonzeros per row in the call to 1991 MatCreateSeqAIJ(). 1992 1993 MUST be called before any calls to MatSetValues(); 1994 1995 @*/ 1996 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1997 { 1998 int ierr,(*f)(Mat,int *); 1999 2000 PetscFunctionBegin; 2001 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2002 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 2003 if (f) { 2004 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2005 } else { 2006 SETERRQ(1,"Wrong type of matrix to set column indices"); 2007 } 2008 PetscFunctionReturn(0); 2009 } 2010 2011 /* ----------------------------------------------------------------------------------------*/ 2012 2013 EXTERN_C_BEGIN 2014 #undef __FUNC__ 2015 #define __FUNC__ /*<a name="MatStoreValues_SeqAIJ"></a>*/"MatStoreValues_SeqAIJ" 2016 int MatStoreValues_SeqAIJ(Mat mat) 2017 { 2018 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2019 int nz = aij->i[aij->m]+aij->indexshift,ierr; 2020 2021 PetscFunctionBegin; 2022 if (aij->nonew != 1) { 2023 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2024 } 2025 2026 /* allocate space for values if not already there */ 2027 if (!aij->saved_values) { 2028 aij->saved_values = (Scalar*)PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 2029 } 2030 2031 /* copy values over */ 2032 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 2033 PetscFunctionReturn(0); 2034 } 2035 EXTERN_C_END 2036 2037 #undef __FUNC__ 2038 #define __FUNC__ /*<a name="MatStoreValues""></a>*/"MatStoreValues" 2039 /*@ 2040 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2041 example, reuse of the linear part of a Jacobian, while recomputing the 2042 nonlinear portion. 2043 2044 Collect on Mat 2045 2046 Input Parameters: 2047 . mat - the matrix (currently on AIJ matrices support this option) 2048 2049 Level: advanced 2050 2051 Common Usage, with SNESSolve(): 2052 $ Create Jacobian matrix 2053 $ Set linear terms into matrix 2054 $ Apply boundary conditions to matrix, at this time matrix must have 2055 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2056 $ boundary conditions again will not change the nonzero structure 2057 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2058 $ ierr = MatStoreValues(mat); 2059 $ Call SNESSetJacobian() with matrix 2060 $ In your Jacobian routine 2061 $ ierr = MatRetrieveValues(mat); 2062 $ Set nonlinear terms in matrix 2063 2064 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2065 $ // build linear portion of Jacobian 2066 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2067 $ ierr = MatStoreValues(mat); 2068 $ loop over nonlinear iterations 2069 $ ierr = MatRetrieveValues(mat); 2070 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2071 $ // call MatAssemblyBegin/End() on matrix 2072 $ Solve linear system with Jacobian 2073 $ endloop 2074 2075 Notes: 2076 Matrix must already be assemblied before calling this routine 2077 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2078 calling this routine. 2079 2080 .seealso: MatRetrieveValues() 2081 2082 @*/ 2083 int MatStoreValues(Mat mat) 2084 { 2085 int ierr,(*f)(Mat); 2086 2087 PetscFunctionBegin; 2088 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2089 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2090 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2091 2092 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f);CHKERRQ(ierr); 2093 if (f) { 2094 ierr = (*f)(mat);CHKERRQ(ierr); 2095 } else { 2096 SETERRQ(1,"Wrong type of matrix to store values"); 2097 } 2098 PetscFunctionReturn(0); 2099 } 2100 2101 EXTERN_C_BEGIN 2102 #undef __FUNC__ 2103 #define __FUNC__ /*<a name="MatRetrieveValues_SeqAIJ"></a>*/"MatRetrieveValues_SeqAIJ" 2104 int MatRetrieveValues_SeqAIJ(Mat mat) 2105 { 2106 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2107 int nz = aij->i[aij->m]+aij->indexshift,ierr; 2108 2109 PetscFunctionBegin; 2110 if (aij->nonew != 1) { 2111 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2112 } 2113 if (!aij->saved_values) { 2114 SETERRQ(1,"Must call MatStoreValues(A);first"); 2115 } 2116 2117 /* copy values over */ 2118 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 2119 PetscFunctionReturn(0); 2120 } 2121 EXTERN_C_END 2122 2123 #undef __FUNC__ 2124 #define __FUNC__ /*<a name="MatRetrieveValues"></a>*/"MatRetrieveValues" 2125 /*@ 2126 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2127 example, reuse of the linear part of a Jacobian, while recomputing the 2128 nonlinear portion. 2129 2130 Collect on Mat 2131 2132 Input Parameters: 2133 . mat - the matrix (currently on AIJ matrices support this option) 2134 2135 Level: advanced 2136 2137 .seealso: MatStoreValues() 2138 2139 @*/ 2140 int MatRetrieveValues(Mat mat) 2141 { 2142 int ierr,(*f)(Mat); 2143 2144 PetscFunctionBegin; 2145 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2146 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2147 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2148 2149 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f);CHKERRQ(ierr); 2150 if (f) { 2151 ierr = (*f)(mat);CHKERRQ(ierr); 2152 } else { 2153 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2154 } 2155 PetscFunctionReturn(0); 2156 } 2157 2158 /* 2159 This allows SeqAIJ matrices to be passed to the matlab engine 2160 */ 2161 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) 2162 #include "engine.h" /* Matlab include file */ 2163 #include "mex.h" /* Matlab include file */ 2164 EXTERN_C_BEGIN 2165 #undef __FUNC__ 2166 #define __FUNC__ /*<a name="MatMatlabEnginePut_SeqAIJ"></a>*/"MatMatlabEnginePut_SeqAIJ" 2167 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2168 { 2169 int ierr,i,*aj,*ai; 2170 Mat B = (Mat)obj; 2171 Scalar *array; 2172 mxArray *mat; 2173 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2174 2175 2176 PetscFunctionBegin; 2177 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2178 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(Scalar));CHKERRQ(ierr); 2179 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2180 ai = mxGetJc(mat); 2181 aj = mxGetIr(mat); 2182 ierr = PetscMemcpy(aj,aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2183 ierr = PetscMemcpy(ai,aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2184 2185 /* Matlab indices start at 0 for sparse (what a surprise) */ 2186 if (aij->indexshift) { 2187 for (i=0; i<B->m+1; i++) { 2188 ai[i]--; 2189 } 2190 for (i=0; i<aij->nz; i++) { 2191 aj[i]--; 2192 } 2193 } 2194 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2195 mxSetName(mat,obj->name); 2196 engPutArray((Engine *)engine,mat); 2197 PetscFunctionReturn(0); 2198 } 2199 EXTERN_C_END 2200 #endif 2201 2202 /* --------------------------------------------------------------------------------*/ 2203 2204 #undef __FUNC__ 2205 #define __FUNC__ /*<a name="MatCreateSeqAIJ"></a>*/"MatCreateSeqAIJ" 2206 /*@C 2207 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2208 (the default parallel PETSc format). For good matrix assembly performance 2209 the user should preallocate the matrix storage by setting the parameter nz 2210 (or the array nnz). By setting these parameters accurately, performance 2211 during matrix assembly can be increased by more than a factor of 50. 2212 2213 Collective on MPI_Comm 2214 2215 Input Parameters: 2216 + comm - MPI communicator, set to PETSC_COMM_SELF 2217 . m - number of rows 2218 . n - number of columns 2219 . nz - number of nonzeros per row (same for all rows) 2220 - nnz - array containing the number of nonzeros in the various rows 2221 (possibly different for each row) or PETSC_NULL 2222 2223 Output Parameter: 2224 . A - the matrix 2225 2226 Notes: 2227 The AIJ format (also called the Yale sparse matrix format or 2228 compressed row storage), is fully compatible with standard Fortran 77 2229 storage. That is, the stored row and column indices can begin at 2230 either one (as in Fortran) or zero. See the users' manual for details. 2231 2232 Specify the preallocated storage with either nz or nnz (not both). 2233 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2234 allocation. For large problems you MUST preallocate memory or you 2235 will get TERRIBLE performance, see the users' manual chapter on matrices. 2236 2237 By default, this format uses inodes (identical nodes) when possible, to 2238 improve numerical efficiency of matrix-vector products and solves. We 2239 search for consecutive rows with the same nonzero structure, thereby 2240 reusing matrix information to achieve increased efficiency. 2241 2242 Options Database Keys: 2243 + -mat_aij_no_inode - Do not use inodes 2244 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2245 - -mat_aij_oneindex - Internally use indexing starting at 1 2246 rather than 0. Note that when calling MatSetValues(), 2247 the user still MUST index entries starting at 0! 2248 2249 Level: intermediate 2250 2251 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2252 2253 @*/ 2254 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2255 { 2256 Mat B; 2257 Mat_SeqAIJ *b; 2258 int i,len,ierr,size; 2259 PetscTruth flg; 2260 2261 PetscFunctionBegin; 2262 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2263 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2264 2265 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 2266 if (nnz) { 2267 for (i=0; i<m; i++) { 2268 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2269 } 2270 } 2271 2272 *A = 0; 2273 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",comm,MatDestroy,MatView); 2274 PLogObjectCreate(B); 2275 B->data = (void*)(b = PetscNew(Mat_SeqAIJ));CHKPTRQ(b); 2276 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2277 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2278 B->factor = 0; 2279 B->lupivotthreshold = 1.0; 2280 B->mapping = 0; 2281 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2282 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2283 b->row = 0; 2284 b->col = 0; 2285 b->icol = 0; 2286 b->indexshift = 0; 2287 b->reallocs = 0; 2288 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2289 if (flg) b->indexshift = -1; 2290 2291 b->m = m; B->m = m; B->M = m; 2292 b->n = n; B->n = n; B->N = n; 2293 2294 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2295 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2296 2297 b->imax = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(b->imax); 2298 if (!nnz) { 2299 if (nz == PETSC_DEFAULT) nz = 10; 2300 else if (nz <= 0) nz = 1; 2301 for (i=0; i<m; i++) b->imax[i] = nz; 2302 nz = nz*m; 2303 } else { 2304 nz = 0; 2305 for (i=0; i<m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2306 } 2307 2308 /* allocate the matrix space */ 2309 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2310 b->a = (Scalar*)PetscMalloc(len);CHKPTRQ(b->a); 2311 b->j = (int*)(b->a + nz); 2312 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2313 b->i = b->j + nz; 2314 b->singlemalloc = PETSC_TRUE; 2315 b->freedata = PETSC_TRUE; 2316 2317 b->i[0] = -b->indexshift; 2318 for (i=1; i<m+1; i++) { 2319 b->i[i] = b->i[i-1] + b->imax[i-1]; 2320 } 2321 2322 /* b->ilen will count nonzeros in each row so far. */ 2323 b->ilen = (int*)PetscMalloc((m+1)*sizeof(int)); 2324 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2325 for (i=0; i<b->m; i++) { b->ilen[i] = 0;} 2326 2327 b->nz = 0; 2328 b->maxnz = nz; 2329 b->sorted = PETSC_FALSE; 2330 b->ignorezeroentries = PETSC_FALSE; 2331 b->roworiented = PETSC_TRUE; 2332 b->nonew = 0; 2333 b->diag = 0; 2334 b->solve_work = 0; 2335 b->spptr = 0; 2336 b->inode.use = PETSC_TRUE; 2337 b->inode.node_count = 0; 2338 b->inode.size = 0; 2339 b->inode.limit = 5; 2340 b->inode.max_limit = 5; 2341 b->saved_values = 0; 2342 B->info.nz_unneeded = (double)b->maxnz; 2343 b->idiag = 0; 2344 b->ssor = 0; 2345 b->keepzeroedrows = PETSC_FALSE; 2346 2347 *A = B; 2348 2349 /* SuperLU is not currently supported through PETSc */ 2350 #if defined(PETSC_HAVE_SUPERLU) 2351 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2352 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2353 #endif 2354 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2355 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2356 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2357 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2358 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2359 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2360 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2361 if (flg) { 2362 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2363 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2364 } 2365 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 2366 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr); } 2367 2368 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2369 "MatSeqAIJSetColumnIndices_SeqAIJ", 2370 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2371 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2372 "MatStoreValues_SeqAIJ", 2373 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2375 "MatRetrieveValues_SeqAIJ", 2376 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2377 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) 2378 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2379 #endif 2380 PetscFunctionReturn(0); 2381 } 2382 2383 #undef __FUNC__ 2384 #define __FUNC__ /*<a name="MatDuplicate_SeqAIJ"></a>*/"MatDuplicate_SeqAIJ" 2385 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2386 { 2387 Mat C; 2388 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2389 int i,len,m = a->m,shift = a->indexshift,ierr; 2390 2391 PetscFunctionBegin; 2392 *B = 0; 2393 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",A->comm,MatDestroy,MatView); 2394 PLogObjectCreate(C); 2395 C->data = (void*)(c = PetscNew(Mat_SeqAIJ));CHKPTRQ(c); 2396 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2397 C->factor = A->factor; 2398 c->row = 0; 2399 c->col = 0; 2400 c->icol = 0; 2401 c->indexshift = shift; 2402 c->keepzeroedrows = a->keepzeroedrows; 2403 C->assembled = PETSC_TRUE; 2404 2405 c->m = C->m = a->m; 2406 c->n = C->n = a->n; 2407 C->M = a->m; 2408 C->N = a->n; 2409 2410 c->imax = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->imax); 2411 c->ilen = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->ilen); 2412 for (i=0; i<m; i++) { 2413 c->imax[i] = a->imax[i]; 2414 c->ilen[i] = a->ilen[i]; 2415 } 2416 2417 /* allocate the matrix space */ 2418 c->singlemalloc = PETSC_TRUE; 2419 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2420 c->a = (Scalar*)PetscMalloc(len);CHKPTRQ(c->a); 2421 c->j = (int*)(c->a + a->i[m] + shift); 2422 c->i = c->j + a->i[m] + shift; 2423 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2424 if (m > 0) { 2425 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2426 if (cpvalues == MAT_COPY_VALUES) { 2427 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2428 } else { 2429 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2430 } 2431 } 2432 2433 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2434 c->sorted = a->sorted; 2435 c->roworiented = a->roworiented; 2436 c->nonew = a->nonew; 2437 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2438 c->saved_values = 0; 2439 c->idiag = 0; 2440 c->ssor = 0; 2441 c->ignorezeroentries = a->ignorezeroentries; 2442 c->freedata = PETSC_TRUE; 2443 2444 if (a->diag) { 2445 c->diag = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->diag); 2446 PLogObjectMemory(C,(m+1)*sizeof(int)); 2447 for (i=0; i<m; i++) { 2448 c->diag[i] = a->diag[i]; 2449 } 2450 } else c->diag = 0; 2451 c->inode.use = a->inode.use; 2452 c->inode.limit = a->inode.limit; 2453 c->inode.max_limit = a->inode.max_limit; 2454 if (a->inode.size){ 2455 c->inode.size = (int*)PetscMalloc((m+1)*sizeof(int));CHKPTRQ(c->inode.size); 2456 c->inode.node_count = a->inode.node_count; 2457 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2458 } else { 2459 c->inode.size = 0; 2460 c->inode.node_count = 0; 2461 } 2462 c->nz = a->nz; 2463 c->maxnz = a->maxnz; 2464 c->solve_work = 0; 2465 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2466 2467 *B = C; 2468 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2469 PetscFunctionReturn(0); 2470 } 2471 2472 #undef __FUNC__ 2473 #define __FUNC__ /*<a name="MatLoad_SeqAIJ"></a>*/"MatLoad_SeqAIJ" 2474 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2475 { 2476 Mat_SeqAIJ *a; 2477 Mat B; 2478 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2479 MPI_Comm comm; 2480 2481 PetscFunctionBegin; 2482 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2483 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2484 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2485 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2486 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2487 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2488 M = header[1]; N = header[2]; nz = header[3]; 2489 2490 if (nz < 0) { 2491 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2492 } 2493 2494 /* read in row lengths */ 2495 rowlengths = (int*)PetscMalloc(M*sizeof(int));CHKPTRQ(rowlengths); 2496 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2497 2498 /* create our matrix */ 2499 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2500 B = *A; 2501 a = (Mat_SeqAIJ*)B->data; 2502 shift = a->indexshift; 2503 2504 /* read in column indices and adjust for Fortran indexing*/ 2505 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2506 if (shift) { 2507 for (i=0; i<nz; i++) { 2508 a->j[i] += 1; 2509 } 2510 } 2511 2512 /* read in nonzero values */ 2513 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2514 2515 /* set matrix "i" values */ 2516 a->i[0] = -shift; 2517 for (i=1; i<= M; i++) { 2518 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2519 a->ilen[i-1] = rowlengths[i-1]; 2520 } 2521 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2522 2523 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2524 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2525 PetscFunctionReturn(0); 2526 } 2527 2528 #undef __FUNC__ 2529 #define __FUNC__ /*<a name="MatEqual_SeqAIJ""></a>*/"MatEqual_SeqAIJ" 2530 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2531 { 2532 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2533 int ierr; 2534 2535 PetscFunctionBegin; 2536 if (B->type != MATSEQAIJ) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2537 2538 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2539 if ((a->m != b->m) || (a->n !=b->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2540 *flg = PETSC_FALSE; 2541 PetscFunctionReturn(0); 2542 } 2543 2544 /* if the a->i are the same */ 2545 ierr = PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2546 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2547 2548 /* if a->j are the same */ 2549 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2550 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2551 2552 /* if a->a are the same */ 2553 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(Scalar),flg);CHKERRQ(ierr); 2554 2555 PetscFunctionReturn(0); 2556 2557 } 2558 2559 #undef __FUNC__ 2560 #define __FUNC__ /*<a name="MatCreateSeqAIJWithArrays"></a>*/"MatCreateSeqAIJWithArrays" 2561 /*@C 2562 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2563 provided by the user. 2564 2565 Coolective on MPI_Comm 2566 2567 Input Parameters: 2568 + comm - must be an MPI communicator of size 1 2569 . m - number of rows 2570 . n - number of columns 2571 . i - row indices 2572 . j - column indices 2573 - a - matrix values 2574 2575 Output Parameter: 2576 . mat - the matrix 2577 2578 Level: intermediate 2579 2580 Notes: 2581 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2582 once the matrix is destroyed 2583 2584 You cannot set new nonzero locations into this matrix, that will generate an error. 2585 2586 The i and j indices can be either 0- or 1 based 2587 2588 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2589 2590 @*/ 2591 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,Scalar *a,Mat *mat) 2592 { 2593 int ierr,ii; 2594 Mat_SeqAIJ *aij; 2595 2596 PetscFunctionBegin; 2597 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2598 aij = (Mat_SeqAIJ*)(*mat)->data; 2599 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2600 2601 if (i[0] == 1) { 2602 aij->indexshift = -1; 2603 } else if (i[0]) { 2604 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2605 } 2606 aij->i = i; 2607 aij->j = j; 2608 aij->a = a; 2609 aij->singlemalloc = PETSC_FALSE; 2610 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2611 aij->freedata = PETSC_FALSE; 2612 2613 for (ii=0; ii<m; ii++) { 2614 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2615 #if defined(PETSC_USE_BOPT_g) 2616 if (i[ii+1] - i[ii] < 0) SETERRQ2(1,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2617 #endif 2618 } 2619 #if defined(PETSC_USE_BOPT_g) 2620 for (ii=0; ii<aij->i[m]; ii++) { 2621 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2622 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2623 } 2624 #endif 2625 2626 /* changes indices to start at 0 */ 2627 if (i[0]) { 2628 aij->indexshift = 0; 2629 for (ii=0; ii<m; ii++) { 2630 i[ii]--; 2631 } 2632 for (ii=0; ii<i[m]; ii++) { 2633 j[ii]--; 2634 } 2635 } 2636 2637 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2638 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2639 PetscFunctionReturn(0); 2640 } 2641 2642 2643 2644