1 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 8 #include "src/inline/spops.h" 9 #include "src/inline/dot.h" 10 #include "petscbt.h" 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 14 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 15 { 16 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 17 PetscErrorCode ierr; 18 PetscInt i,ishift; 19 20 PetscFunctionBegin; 21 *m = A->m; 22 if (!ia) PetscFunctionReturn(0); 23 ishift = 0; 24 if (symmetric && !A->structurally_symmetric) { 25 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 26 } else if (oshift == 1) { 27 PetscInt nz = a->i[A->m]; 28 /* malloc space and add 1 to i and j indices */ 29 ierr = PetscMalloc((A->m+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 30 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 31 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 32 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 33 } else { 34 *ia = a->i; *ja = a->j; 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 41 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 42 { 43 PetscErrorCode ierr; 44 45 PetscFunctionBegin; 46 if (!ia) PetscFunctionReturn(0); 47 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 48 ierr = PetscFree(*ia);CHKERRQ(ierr); 49 ierr = PetscFree(*ja);CHKERRQ(ierr); 50 } 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 56 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 57 { 58 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 59 PetscErrorCode ierr; 60 PetscInt i,*collengths,*cia,*cja,n = A->n,m = A->m; 61 PetscInt nz = a->i[m],row,*jj,mr,col; 62 63 PetscFunctionBegin; 64 *nn = A->n; 65 if (!ia) PetscFunctionReturn(0); 66 if (symmetric) { 67 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 68 } else { 69 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 70 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 71 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 72 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 73 jj = a->j; 74 for (i=0; i<nz; i++) { 75 collengths[jj[i]]++; 76 } 77 cia[0] = oshift; 78 for (i=0; i<n; i++) { 79 cia[i+1] = cia[i] + collengths[i]; 80 } 81 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 82 jj = a->j; 83 for (row=0; row<m; row++) { 84 mr = a->i[row+1] - a->i[row]; 85 for (i=0; i<mr; i++) { 86 col = *jj++; 87 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 88 } 89 } 90 ierr = PetscFree(collengths);CHKERRQ(ierr); 91 *ia = cia; *ja = cja; 92 } 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 98 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 99 { 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 if (!ia) PetscFunctionReturn(0); 104 105 ierr = PetscFree(*ia);CHKERRQ(ierr); 106 ierr = PetscFree(*ja);CHKERRQ(ierr); 107 108 PetscFunctionReturn(0); 109 } 110 111 #define CHUNKSIZE 15 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatSetValues_SeqAIJ" 115 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 116 { 117 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 118 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted = a->sorted; 119 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 120 PetscErrorCode ierr; 121 PetscInt *aj = a->j,nonew = a->nonew; 122 PetscScalar *ap,value,*aa = a->a; 123 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 124 PetscTruth roworiented = a->roworiented; 125 126 PetscFunctionBegin; 127 for (k=0; k<m; k++) { /* loop over added rows */ 128 row = im[k]; 129 if (row < 0) continue; 130 #if defined(PETSC_USE_BOPT_g) 131 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 132 #endif 133 rp = aj + ai[row]; ap = aa + ai[row]; 134 rmax = imax[row]; nrow = ailen[row]; 135 low = 0; 136 for (l=0; l<n; l++) { /* loop over added columns */ 137 if (in[l] < 0) continue; 138 #if defined(PETSC_USE_BOPT_g) 139 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 140 #endif 141 col = in[l]; 142 if (roworiented) { 143 value = v[l + k*n]; 144 } else { 145 value = v[k + l*m]; 146 } 147 if (value == 0.0 && ignorezeroentries) continue; 148 149 if (!sorted) low = 0; high = nrow; 150 while (high-low > 5) { 151 t = (low+high)/2; 152 if (rp[t] > col) high = t; 153 else low = t; 154 } 155 for (i=low; i<high; i++) { 156 if (rp[i] > col) break; 157 if (rp[i] == col) { 158 if (is == ADD_VALUES) ap[i] += value; 159 else ap[i] = value; 160 goto noinsert; 161 } 162 } 163 if (nonew == 1) goto noinsert; 164 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 165 if (nrow >= rmax) { 166 /* there is no extra room in row, therefore enlarge */ 167 PetscInt new_nz = ai[A->m] + CHUNKSIZE,*new_i,*new_j; 168 size_t len; 169 PetscScalar *new_a; 170 171 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix requiring new malloc()",row,col); 172 173 /* malloc new storage space */ 174 len = ((size_t) new_nz)*(sizeof(PetscInt)+sizeof(PetscScalar))+(A->m+1)*sizeof(PetscInt); 175 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 176 new_j = (PetscInt*)(new_a + new_nz); 177 new_i = new_j + new_nz; 178 179 /* copy over old data into new slots */ 180 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 181 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 182 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(PetscInt));CHKERRQ(ierr); 183 len = (((size_t) new_nz) - CHUNKSIZE - ai[row] - nrow ); 184 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(PetscInt));CHKERRQ(ierr); 185 ierr = PetscMemcpy(new_a,aa,(((size_t) ai[row])+nrow)*sizeof(PetscScalar));CHKERRQ(ierr); 186 ierr = PetscMemcpy(new_a+ai[row]+nrow+CHUNKSIZE,aa+ai[row]+nrow,len*sizeof(PetscScalar));CHKERRQ(ierr); 187 /* free up old matrix storage */ 188 ierr = PetscFree(a->a);CHKERRQ(ierr); 189 if (!a->singlemalloc) { 190 ierr = PetscFree(a->i);CHKERRQ(ierr); 191 ierr = PetscFree(a->j);CHKERRQ(ierr); 192 } 193 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 194 a->singlemalloc = PETSC_TRUE; 195 196 rp = aj + ai[row]; ap = aa + ai[row] ; 197 rmax = imax[row] = imax[row] + CHUNKSIZE; 198 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(PetscInt) + sizeof(PetscScalar))); 199 a->maxnz += CHUNKSIZE; 200 a->reallocs++; 201 } 202 N = nrow++ - 1; a->nz++; 203 /* shift up all the later entries in this row */ 204 for (ii=N; ii>=i; ii--) { 205 rp[ii+1] = rp[ii]; 206 ap[ii+1] = ap[ii]; 207 } 208 rp[i] = col; 209 ap[i] = value; 210 noinsert:; 211 low = i + 1; 212 } 213 ailen[row] = nrow; 214 } 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "MatGetValues_SeqAIJ" 220 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 221 { 222 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 223 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 224 PetscInt *ai = a->i,*ailen = a->ilen; 225 PetscScalar *ap,*aa = a->a,zero = 0.0; 226 227 PetscFunctionBegin; 228 for (k=0; k<m; k++) { /* loop over rows */ 229 row = im[k]; 230 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 231 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 232 rp = aj + ai[row]; ap = aa + ai[row]; 233 nrow = ailen[row]; 234 for (l=0; l<n; l++) { /* loop over columns */ 235 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 236 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 237 col = in[l] ; 238 high = nrow; low = 0; /* assume unsorted */ 239 while (high-low > 5) { 240 t = (low+high)/2; 241 if (rp[t] > col) high = t; 242 else low = t; 243 } 244 for (i=low; i<high; i++) { 245 if (rp[i] > col) break; 246 if (rp[i] == col) { 247 *v++ = ap[i]; 248 goto finished; 249 } 250 } 251 *v++ = zero; 252 finished:; 253 } 254 } 255 PetscFunctionReturn(0); 256 } 257 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "MatView_SeqAIJ_Binary" 261 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 262 { 263 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 264 PetscErrorCode ierr; 265 PetscInt i,*col_lens; 266 int fd; 267 268 PetscFunctionBegin; 269 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 270 ierr = PetscMalloc((4+A->m)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 271 col_lens[0] = MAT_FILE_COOKIE; 272 col_lens[1] = A->m; 273 col_lens[2] = A->n; 274 col_lens[3] = a->nz; 275 276 /* store lengths of each row and write (including header) to file */ 277 for (i=0; i<A->m; i++) { 278 col_lens[4+i] = a->i[i+1] - a->i[i]; 279 } 280 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 281 ierr = PetscFree(col_lens);CHKERRQ(ierr); 282 283 /* store column indices (zero start index) */ 284 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 285 286 /* store nonzero values */ 287 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 295 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 296 { 297 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 298 PetscErrorCode ierr; 299 PetscInt i,j,m = A->m,shift=0; 300 char *name; 301 PetscViewerFormat format; 302 303 PetscFunctionBegin; 304 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 305 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 306 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_INFO) { 307 if (a->inode.size) { 308 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %D nodes, limit used is %D\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 309 } else { 310 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 311 } 312 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 313 PetscInt nofinalvalue = 0; 314 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 315 nofinalvalue = 1; 316 } 317 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 318 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->n);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 322 323 for (i=0; i<m; i++) { 324 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 325 #if defined(PETSC_USE_COMPLEX) 326 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 327 #else 328 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 329 #endif 330 } 331 } 332 if (nofinalvalue) { 333 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 334 } 335 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 337 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 338 PetscFunctionReturn(0); 339 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 340 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 341 for (i=0; i<m; i++) { 342 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 343 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 344 #if defined(PETSC_USE_COMPLEX) 345 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 346 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 347 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 348 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 349 } else if (PetscRealPart(a->a[j]) != 0.0) { 350 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 351 } 352 #else 353 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 354 #endif 355 } 356 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 357 } 358 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 359 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 360 PetscInt nzd=0,fshift=1,*sptr; 361 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 362 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 363 for (i=0; i<m; i++) { 364 sptr[i] = nzd+1; 365 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 366 if (a->j[j] >= i) { 367 #if defined(PETSC_USE_COMPLEX) 368 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 369 #else 370 if (a->a[j] != 0.0) nzd++; 371 #endif 372 } 373 } 374 } 375 sptr[m] = nzd+1; 376 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 377 for (i=0; i<m+1; i+=6) { 378 if (i+4<m) {ierr = PetscViewerASCIIPrintf(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);} 379 else if (i+3<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]);CHKERRQ(ierr);} 380 else if (i+2<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]);CHKERRQ(ierr);} 381 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 382 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 383 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 384 } 385 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 386 ierr = PetscFree(sptr);CHKERRQ(ierr); 387 for (i=0; i<m; i++) { 388 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 389 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 390 } 391 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 392 } 393 ierr = PetscViewerASCIIPrintf(viewer,"\n");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) { 397 #if defined(PETSC_USE_COMPLEX) 398 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 399 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 400 } 401 #else 402 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 403 #endif 404 } 405 } 406 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 407 } 408 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 409 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 410 PetscInt cnt = 0,jcnt; 411 PetscScalar value; 412 413 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 414 for (i=0; i<m; i++) { 415 jcnt = 0; 416 for (j=0; j<A->n; j++) { 417 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 418 value = a->a[cnt++]; 419 jcnt++; 420 } else { 421 value = 0.0; 422 } 423 #if defined(PETSC_USE_COMPLEX) 424 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 425 #else 426 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 427 #endif 428 } 429 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 430 } 431 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 432 } else { 433 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 434 for (i=0; i<m; i++) { 435 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 436 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 437 #if defined(PETSC_USE_COMPLEX) 438 if (PetscImaginaryPart(a->a[j]) > 0.0) { 439 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 440 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 441 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 442 } else { 443 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 444 } 445 #else 446 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 447 #endif 448 } 449 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 450 } 451 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 452 } 453 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 454 PetscFunctionReturn(0); 455 } 456 457 #undef __FUNCT__ 458 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 459 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 460 { 461 Mat A = (Mat) Aa; 462 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 463 PetscErrorCode ierr; 464 PetscInt i,j,m = A->m,color; 465 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 466 PetscViewer viewer; 467 PetscViewerFormat format; 468 469 PetscFunctionBegin; 470 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 471 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 472 473 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 474 /* loop over matrix elements drawing boxes */ 475 476 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 477 /* Blue for negative, Cyan for zero and Red for positive */ 478 color = PETSC_DRAW_BLUE; 479 for (i=0; i<m; i++) { 480 y_l = m - i - 1.0; y_r = y_l + 1.0; 481 for (j=a->i[i]; j<a->i[i+1]; j++) { 482 x_l = a->j[j] ; x_r = x_l + 1.0; 483 #if defined(PETSC_USE_COMPLEX) 484 if (PetscRealPart(a->a[j]) >= 0.) continue; 485 #else 486 if (a->a[j] >= 0.) continue; 487 #endif 488 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 489 } 490 } 491 color = PETSC_DRAW_CYAN; 492 for (i=0; i<m; i++) { 493 y_l = m - i - 1.0; y_r = y_l + 1.0; 494 for (j=a->i[i]; j<a->i[i+1]; j++) { 495 x_l = a->j[j]; x_r = x_l + 1.0; 496 if (a->a[j] != 0.) continue; 497 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 498 } 499 } 500 color = PETSC_DRAW_RED; 501 for (i=0; i<m; i++) { 502 y_l = m - i - 1.0; y_r = y_l + 1.0; 503 for (j=a->i[i]; j<a->i[i+1]; j++) { 504 x_l = a->j[j]; x_r = x_l + 1.0; 505 #if defined(PETSC_USE_COMPLEX) 506 if (PetscRealPart(a->a[j]) <= 0.) continue; 507 #else 508 if (a->a[j] <= 0.) continue; 509 #endif 510 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 511 } 512 } 513 } else { 514 /* use contour shading to indicate magnitude of values */ 515 /* first determine max of all nonzero values */ 516 PetscInt nz = a->nz,count; 517 PetscDraw popup; 518 PetscReal scale; 519 520 for (i=0; i<nz; i++) { 521 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 522 } 523 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 524 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 525 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 526 count = 0; 527 for (i=0; i<m; i++) { 528 y_l = m - i - 1.0; y_r = y_l + 1.0; 529 for (j=a->i[i]; j<a->i[i+1]; j++) { 530 x_l = a->j[j]; x_r = x_l + 1.0; 531 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 532 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 533 count++; 534 } 535 } 536 } 537 PetscFunctionReturn(0); 538 } 539 540 #undef __FUNCT__ 541 #define __FUNCT__ "MatView_SeqAIJ_Draw" 542 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 543 { 544 PetscErrorCode ierr; 545 PetscDraw draw; 546 PetscReal xr,yr,xl,yl,h,w; 547 PetscTruth isnull; 548 549 PetscFunctionBegin; 550 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 551 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 552 if (isnull) PetscFunctionReturn(0); 553 554 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 555 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 556 xr += w; yr += h; xl = -w; yl = -h; 557 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 558 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 559 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 560 PetscFunctionReturn(0); 561 } 562 563 #undef __FUNCT__ 564 #define __FUNCT__ "MatView_SeqAIJ" 565 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 566 { 567 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 568 PetscErrorCode ierr; 569 PetscTruth issocket,iascii,isbinary,isdraw; 570 571 PetscFunctionBegin; 572 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 573 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 574 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 575 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 576 if (issocket) { 577 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 578 } else if (iascii) { 579 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 580 } else if (isbinary) { 581 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 582 } else if (isdraw) { 583 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 584 } else { 585 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 586 } 587 PetscFunctionReturn(0); 588 } 589 590 #undef __FUNCT__ 591 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 592 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 595 PetscErrorCode ierr; 596 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 597 PetscInt m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 598 PetscScalar *aa = a->a,*ap; 599 600 PetscFunctionBegin; 601 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 602 603 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 604 for (i=1; i<m; i++) { 605 /* move each row back by the amount of empty slots (fshift) before it*/ 606 fshift += imax[i-1] - ailen[i-1]; 607 rmax = PetscMax(rmax,ailen[i]); 608 if (fshift) { 609 ip = aj + ai[i] ; 610 ap = aa + ai[i] ; 611 N = ailen[i]; 612 for (j=0; j<N; j++) { 613 ip[j-fshift] = ip[j]; 614 ap[j-fshift] = ap[j]; 615 } 616 } 617 ai[i] = ai[i-1] + ailen[i-1]; 618 } 619 if (m) { 620 fshift += imax[m-1] - ailen[m-1]; 621 ai[m] = ai[m-1] + ailen[m-1]; 622 } 623 /* reset ilen and imax for each row */ 624 for (i=0; i<m; i++) { 625 ailen[i] = imax[i] = ai[i+1] - ai[i]; 626 } 627 a->nz = ai[m]; 628 629 /* diagonals may have moved, so kill the diagonal pointers */ 630 if (fshift && a->diag) { 631 ierr = PetscFree(a->diag);CHKERRQ(ierr); 632 PetscLogObjectMemory(A,-(m+1)*sizeof(PetscInt)); 633 a->diag = 0; 634 } 635 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->n,fshift,a->nz); 636 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %D\n",a->reallocs); 637 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Maximum nonzeros in any row is %D\n",rmax); 638 a->reallocs = 0; 639 A->info.nz_unneeded = (double)fshift; 640 a->rmax = rmax; 641 642 /* check out for identical nodes. If found, use inode functions */ 643 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 644 645 /* check out for zero rows. If found, use CompressedRow functions */ 646 ierr = Mat_AIJ_CheckCompressedRow(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 647 648 PetscFunctionReturn(0); 649 } 650 651 #undef __FUNCT__ 652 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 653 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 654 { 655 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 660 PetscFunctionReturn(0); 661 } 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "MatDestroy_SeqAIJ" 665 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 666 { 667 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 #if defined(PETSC_USE_LOG) 672 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->m,A->n,a->nz); 673 #endif 674 if (a->freedata) { 675 ierr = PetscFree(a->a);CHKERRQ(ierr); 676 if (!a->singlemalloc) { 677 ierr = PetscFree(a->i);CHKERRQ(ierr); 678 ierr = PetscFree(a->j);CHKERRQ(ierr); 679 } 680 } 681 if (a->row) { 682 ierr = ISDestroy(a->row);CHKERRQ(ierr); 683 } 684 if (a->col) { 685 ierr = ISDestroy(a->col);CHKERRQ(ierr); 686 } 687 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 688 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 689 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 690 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 691 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 692 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 693 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 694 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 695 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 696 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 697 if (a->compressedrow.use){ierr = PetscFree(a->compressedrow.i);} 698 699 ierr = PetscFree(a);CHKERRQ(ierr); 700 701 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 702 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 703 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 704 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 706 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 707 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 709 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatAdjustForInodes_C","",PETSC_NULL);CHKERRQ(ierr); 710 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJGetInodeSizes_C","",PETSC_NULL);CHKERRQ(ierr); 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNCT__ 715 #define __FUNCT__ "MatCompress_SeqAIJ" 716 PetscErrorCode MatCompress_SeqAIJ(Mat A) 717 { 718 PetscFunctionBegin; 719 PetscFunctionReturn(0); 720 } 721 722 #undef __FUNCT__ 723 #define __FUNCT__ "MatSetOption_SeqAIJ" 724 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op) 725 { 726 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 727 728 PetscFunctionBegin; 729 switch (op) { 730 case MAT_ROW_ORIENTED: 731 a->roworiented = PETSC_TRUE; 732 break; 733 case MAT_KEEP_ZEROED_ROWS: 734 a->keepzeroedrows = PETSC_TRUE; 735 break; 736 case MAT_COLUMN_ORIENTED: 737 a->roworiented = PETSC_FALSE; 738 break; 739 case MAT_COLUMNS_SORTED: 740 a->sorted = PETSC_TRUE; 741 break; 742 case MAT_COLUMNS_UNSORTED: 743 a->sorted = PETSC_FALSE; 744 break; 745 case MAT_NO_NEW_NONZERO_LOCATIONS: 746 a->nonew = 1; 747 break; 748 case MAT_NEW_NONZERO_LOCATION_ERR: 749 a->nonew = -1; 750 break; 751 case MAT_NEW_NONZERO_ALLOCATION_ERR: 752 a->nonew = -2; 753 break; 754 case MAT_YES_NEW_NONZERO_LOCATIONS: 755 a->nonew = 0; 756 break; 757 case MAT_IGNORE_ZERO_ENTRIES: 758 a->ignorezeroentries = PETSC_TRUE; 759 break; 760 case MAT_USE_INODES: 761 a->inode.use = PETSC_TRUE; 762 break; 763 case MAT_DO_NOT_USE_INODES: 764 a->inode.use = PETSC_FALSE; 765 break; 766 case MAT_USE_COMPRESSEDROW: 767 a->compressedrow.use = PETSC_TRUE; 768 break; 769 case MAT_DO_NOT_USE_COMPRESSEDROW: 770 a->compressedrow.use = PETSC_FALSE; 771 break; 772 case MAT_ROWS_SORTED: 773 case MAT_ROWS_UNSORTED: 774 case MAT_YES_NEW_DIAGONALS: 775 case MAT_IGNORE_OFF_PROC_ENTRIES: 776 case MAT_USE_HASH_TABLE: 777 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 778 break; 779 case MAT_NO_NEW_DIAGONALS: 780 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 781 case MAT_INODE_LIMIT_1: 782 a->inode.limit = 1; 783 break; 784 case MAT_INODE_LIMIT_2: 785 a->inode.limit = 2; 786 break; 787 case MAT_INODE_LIMIT_3: 788 a->inode.limit = 3; 789 break; 790 case MAT_INODE_LIMIT_4: 791 a->inode.limit = 4; 792 break; 793 case MAT_INODE_LIMIT_5: 794 a->inode.limit = 5; 795 break; 796 case MAT_SYMMETRIC: 797 case MAT_STRUCTURALLY_SYMMETRIC: 798 case MAT_NOT_SYMMETRIC: 799 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 800 case MAT_HERMITIAN: 801 case MAT_NOT_HERMITIAN: 802 case MAT_SYMMETRY_ETERNAL: 803 case MAT_NOT_SYMMETRY_ETERNAL: 804 break; 805 default: 806 SETERRQ(PETSC_ERR_SUP,"unknown option"); 807 } 808 PetscFunctionReturn(0); 809 } 810 811 #undef __FUNCT__ 812 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 813 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 814 { 815 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 816 PetscErrorCode ierr; 817 PetscInt i,j,n; 818 PetscScalar *x,zero = 0.0; 819 820 PetscFunctionBegin; 821 ierr = VecSet(&zero,v);CHKERRQ(ierr); 822 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 823 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 824 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 825 for (i=0; i<A->m; i++) { 826 for (j=a->i[i]; j<a->i[i+1]; j++) { 827 if (a->j[j] == i) { 828 x[i] = a->a[j]; 829 break; 830 } 831 } 832 } 833 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 834 PetscFunctionReturn(0); 835 } 836 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 840 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 841 { 842 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 843 PetscScalar *x,*y; 844 PetscErrorCode ierr; 845 PetscInt m = A->m; 846 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 847 PetscScalar *v,alpha; 848 PetscInt n,i,*idx; 849 #endif 850 851 PetscFunctionBegin; 852 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 853 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 854 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 855 856 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 857 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 858 #else 859 for (i=0; i<m; i++) { 860 idx = a->j + a->i[i] ; 861 v = a->a + a->i[i] ; 862 n = a->i[i+1] - a->i[i]; 863 alpha = x[i]; 864 while (n-->0) {y[*idx++] += alpha * *v++;} 865 } 866 #endif 867 PetscLogFlops(2*a->nz); 868 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 869 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 870 PetscFunctionReturn(0); 871 } 872 873 #undef __FUNCT__ 874 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 875 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 876 { 877 PetscScalar zero = 0.0; 878 PetscErrorCode ierr; 879 880 PetscFunctionBegin; 881 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 882 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 887 #undef __FUNCT__ 888 #define __FUNCT__ "MatMult_SeqAIJ" 889 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 890 { 891 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 892 PetscScalar *x,*y,*v; 893 PetscErrorCode ierr; 894 PetscInt m = A->m,*idx,*ii; 895 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 896 PetscInt n,i,jrow,j; 897 PetscScalar sum; 898 #endif 899 900 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 901 #pragma disjoint(*x,*y,*v) 902 #endif 903 904 PetscFunctionBegin; 905 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 906 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 907 idx = a->j; 908 v = a->a; 909 ii = a->i; 910 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 911 fortranmultaij_(&m,x,ii,idx,v,y); 912 #else 913 for (i=0; i<m; i++) { 914 jrow = ii[i]; 915 n = ii[i+1] - jrow; 916 sum = 0.0; 917 for (j=0; j<n; j++) { 918 sum += v[jrow]*x[idx[jrow]]; jrow++; 919 } 920 y[i] = sum; 921 } 922 #endif 923 PetscLogFlops(2*a->nz - m); 924 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 925 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 926 PetscFunctionReturn(0); 927 } 928 929 #undef __FUNCT__ 930 #define __FUNCT__ "MatMultAdd_SeqAIJ" 931 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 932 { 933 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 934 PetscScalar *x,*y,*z,*v; 935 PetscErrorCode ierr; 936 PetscInt m = A->m,*idx,*ii; 937 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 938 PetscInt n,i,jrow,j; 939 PetscScalar sum; 940 #endif 941 942 PetscFunctionBegin; 943 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 944 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 945 if (zz != yy) { 946 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 947 } else { 948 z = y; 949 } 950 951 idx = a->j; 952 v = a->a; 953 ii = a->i; 954 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 955 fortranmultaddaij_(&m,x,ii,idx,v,y,z); 956 #else 957 for (i=0; i<m; i++) { 958 jrow = ii[i]; 959 n = ii[i+1] - jrow; 960 sum = y[i]; 961 for (j=0; j<n; j++) { 962 sum += v[jrow]*x[idx[jrow]]; jrow++; 963 } 964 z[i] = sum; 965 } 966 #endif 967 PetscLogFlops(2*a->nz); 968 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 969 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 970 if (zz != yy) { 971 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 972 } 973 PetscFunctionReturn(0); 974 } 975 976 /* 977 Adds diagonal pointers to sparse matrix structure. 978 */ 979 #undef __FUNCT__ 980 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 981 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 982 { 983 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 984 PetscErrorCode ierr; 985 PetscInt i,j,*diag,m = A->m; 986 987 PetscFunctionBegin; 988 if (a->diag) PetscFunctionReturn(0); 989 990 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&diag);CHKERRQ(ierr); 991 PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt)); 992 for (i=0; i<A->m; i++) { 993 diag[i] = a->i[i+1]; 994 for (j=a->i[i]; j<a->i[i+1]; j++) { 995 if (a->j[j] == i) { 996 diag[i] = j; 997 break; 998 } 999 } 1000 } 1001 a->diag = diag; 1002 PetscFunctionReturn(0); 1003 } 1004 1005 /* 1006 Checks for missing diagonals 1007 */ 1008 #undef __FUNCT__ 1009 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1010 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A) 1011 { 1012 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1013 PetscErrorCode ierr; 1014 PetscInt *diag,*jj = a->j,i; 1015 1016 PetscFunctionBegin; 1017 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1018 diag = a->diag; 1019 for (i=0; i<A->m; i++) { 1020 if (jj[diag[i]] != i) { 1021 SETERRQ1(PETSC_ERR_PLIB,"Matrix is missing diagonal number %D",i); 1022 } 1023 } 1024 PetscFunctionReturn(0); 1025 } 1026 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "MatRelax_SeqAIJ" 1029 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1030 { 1031 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1032 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1033 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1034 PetscErrorCode ierr; 1035 PetscInt n = A->n,m = A->m,i; 1036 const PetscInt *idx,*diag; 1037 1038 PetscFunctionBegin; 1039 its = its*lits; 1040 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1041 1042 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1043 diag = a->diag; 1044 if (!a->idiag) { 1045 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1046 a->ssor = a->idiag + m; 1047 mdiag = a->ssor + m; 1048 1049 v = a->a; 1050 1051 /* this is wrong when fshift omega changes each iteration */ 1052 if (omega == 1.0 && !fshift) { 1053 for (i=0; i<m; i++) { 1054 mdiag[i] = v[diag[i]]; 1055 a->idiag[i] = 1.0/v[diag[i]]; 1056 } 1057 PetscLogFlops(m); 1058 } else { 1059 for (i=0; i<m; i++) { 1060 mdiag[i] = v[diag[i]]; 1061 a->idiag[i] = omega/(fshift + v[diag[i]]); 1062 } 1063 PetscLogFlops(2*m); 1064 } 1065 } 1066 t = a->ssor; 1067 idiag = a->idiag; 1068 mdiag = a->idiag + 2*m; 1069 1070 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1071 if (xx != bb) { 1072 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1073 } else { 1074 b = x; 1075 } 1076 1077 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1078 xs = x; 1079 if (flag == SOR_APPLY_UPPER) { 1080 /* apply (U + D/omega) to the vector */ 1081 bs = b; 1082 for (i=0; i<m; i++) { 1083 d = fshift + a->a[diag[i]]; 1084 n = a->i[i+1] - diag[i] - 1; 1085 idx = a->j + diag[i] + 1; 1086 v = a->a + diag[i] + 1; 1087 sum = b[i]*d/omega; 1088 SPARSEDENSEDOT(sum,bs,v,idx,n); 1089 x[i] = sum; 1090 } 1091 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1092 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1093 PetscLogFlops(a->nz); 1094 PetscFunctionReturn(0); 1095 } 1096 1097 1098 /* Let A = L + U + D; where L is lower trianglar, 1099 U is upper triangular, E is diagonal; This routine applies 1100 1101 (L + E)^{-1} A (U + E)^{-1} 1102 1103 to a vector efficiently using Eisenstat's trick. This is for 1104 the case of SSOR preconditioner, so E is D/omega where omega 1105 is the relaxation factor. 1106 */ 1107 1108 if (flag == SOR_APPLY_LOWER) { 1109 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1110 } else if (flag & SOR_EISENSTAT) { 1111 /* Let A = L + U + D; where L is lower trianglar, 1112 U is upper triangular, E is diagonal; This routine applies 1113 1114 (L + E)^{-1} A (U + E)^{-1} 1115 1116 to a vector efficiently using Eisenstat's trick. This is for 1117 the case of SSOR preconditioner, so E is D/omega where omega 1118 is the relaxation factor. 1119 */ 1120 scale = (2.0/omega) - 1.0; 1121 1122 /* x = (E + U)^{-1} b */ 1123 for (i=m-1; i>=0; i--) { 1124 n = a->i[i+1] - diag[i] - 1; 1125 idx = a->j + diag[i] + 1; 1126 v = a->a + diag[i] + 1; 1127 sum = b[i]; 1128 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1129 x[i] = sum*idiag[i]; 1130 } 1131 1132 /* t = b - (2*E - D)x */ 1133 v = a->a; 1134 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1135 1136 /* t = (E + L)^{-1}t */ 1137 ts = t; 1138 diag = a->diag; 1139 for (i=0; i<m; i++) { 1140 n = diag[i] - a->i[i]; 1141 idx = a->j + a->i[i]; 1142 v = a->a + a->i[i]; 1143 sum = t[i]; 1144 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1145 t[i] = sum*idiag[i]; 1146 /* x = x + t */ 1147 x[i] += t[i]; 1148 } 1149 1150 PetscLogFlops(6*m-1 + 2*a->nz); 1151 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1152 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1153 PetscFunctionReturn(0); 1154 } 1155 if (flag & SOR_ZERO_INITIAL_GUESS) { 1156 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1157 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1158 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1159 #else 1160 for (i=0; i<m; i++) { 1161 n = diag[i] - a->i[i]; 1162 idx = a->j + a->i[i]; 1163 v = a->a + a->i[i]; 1164 sum = b[i]; 1165 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1166 x[i] = sum*idiag[i]; 1167 } 1168 #endif 1169 xb = x; 1170 PetscLogFlops(a->nz); 1171 } else xb = b; 1172 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1173 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1174 for (i=0; i<m; i++) { 1175 x[i] *= mdiag[i]; 1176 } 1177 PetscLogFlops(m); 1178 } 1179 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1180 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1181 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1182 #else 1183 for (i=m-1; i>=0; i--) { 1184 n = a->i[i+1] - diag[i] - 1; 1185 idx = a->j + diag[i] + 1; 1186 v = a->a + diag[i] + 1; 1187 sum = xb[i]; 1188 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1189 x[i] = sum*idiag[i]; 1190 } 1191 #endif 1192 PetscLogFlops(a->nz); 1193 } 1194 its--; 1195 } 1196 while (its--) { 1197 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1198 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1199 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1200 #else 1201 for (i=0; i<m; i++) { 1202 d = fshift + a->a[diag[i]]; 1203 n = a->i[i+1] - a->i[i]; 1204 idx = a->j + a->i[i]; 1205 v = a->a + a->i[i]; 1206 sum = b[i]; 1207 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1208 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1209 } 1210 #endif 1211 PetscLogFlops(a->nz); 1212 } 1213 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1214 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1215 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1216 #else 1217 for (i=m-1; i>=0; i--) { 1218 d = fshift + a->a[diag[i]]; 1219 n = a->i[i+1] - a->i[i]; 1220 idx = a->j + a->i[i]; 1221 v = a->a + a->i[i]; 1222 sum = b[i]; 1223 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1224 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1225 } 1226 #endif 1227 PetscLogFlops(a->nz); 1228 } 1229 } 1230 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1231 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1232 PetscFunctionReturn(0); 1233 } 1234 1235 #undef __FUNCT__ 1236 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1237 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1238 { 1239 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1240 1241 PetscFunctionBegin; 1242 info->rows_global = (double)A->m; 1243 info->columns_global = (double)A->n; 1244 info->rows_local = (double)A->m; 1245 info->columns_local = (double)A->n; 1246 info->block_size = 1.0; 1247 info->nz_allocated = (double)a->maxnz; 1248 info->nz_used = (double)a->nz; 1249 info->nz_unneeded = (double)(a->maxnz - a->nz); 1250 info->assemblies = (double)A->num_ass; 1251 info->mallocs = (double)a->reallocs; 1252 info->memory = A->mem; 1253 if (A->factor) { 1254 info->fill_ratio_given = A->info.fill_ratio_given; 1255 info->fill_ratio_needed = A->info.fill_ratio_needed; 1256 info->factor_mallocs = A->info.factor_mallocs; 1257 } else { 1258 info->fill_ratio_given = 0; 1259 info->fill_ratio_needed = 0; 1260 info->factor_mallocs = 0; 1261 } 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1267 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1268 { 1269 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1270 PetscErrorCode ierr; 1271 PetscInt i,N,*rows,m = A->m - 1; 1272 1273 PetscFunctionBegin; 1274 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1275 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1276 if (a->keepzeroedrows) { 1277 for (i=0; i<N; i++) { 1278 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1279 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1280 } 1281 if (diag) { 1282 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1283 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1284 for (i=0; i<N; i++) { 1285 a->a[a->diag[rows[i]]] = *diag; 1286 } 1287 } 1288 } else { 1289 if (diag) { 1290 for (i=0; i<N; i++) { 1291 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1292 if (a->ilen[rows[i]] > 0) { 1293 a->ilen[rows[i]] = 1; 1294 a->a[a->i[rows[i]]] = *diag; 1295 a->j[a->i[rows[i]]] = rows[i]; 1296 } else { /* in case row was completely empty */ 1297 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1298 } 1299 } 1300 } else { 1301 for (i=0; i<N; i++) { 1302 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1303 a->ilen[rows[i]] = 0; 1304 } 1305 } 1306 } 1307 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1308 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 #undef __FUNCT__ 1313 #define __FUNCT__ "MatGetRow_SeqAIJ" 1314 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1315 { 1316 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1317 PetscInt *itmp; 1318 1319 PetscFunctionBegin; 1320 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1321 1322 *nz = a->i[row+1] - a->i[row]; 1323 if (v) *v = a->a + a->i[row]; 1324 if (idx) { 1325 itmp = a->j + a->i[row]; 1326 if (*nz) { 1327 *idx = itmp; 1328 } 1329 else *idx = 0; 1330 } 1331 PetscFunctionReturn(0); 1332 } 1333 1334 /* remove this function? */ 1335 #undef __FUNCT__ 1336 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1337 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1338 { 1339 PetscFunctionBegin; 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "MatNorm_SeqAIJ" 1345 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1346 { 1347 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1348 PetscScalar *v = a->a; 1349 PetscReal sum = 0.0; 1350 PetscErrorCode ierr; 1351 PetscInt i,j; 1352 1353 PetscFunctionBegin; 1354 if (type == NORM_FROBENIUS) { 1355 for (i=0; i<a->nz; i++) { 1356 #if defined(PETSC_USE_COMPLEX) 1357 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1358 #else 1359 sum += (*v)*(*v); v++; 1360 #endif 1361 } 1362 *nrm = sqrt(sum); 1363 } else if (type == NORM_1) { 1364 PetscReal *tmp; 1365 PetscInt *jj = a->j; 1366 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1367 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1368 *nrm = 0.0; 1369 for (j=0; j<a->nz; j++) { 1370 tmp[*jj++] += PetscAbsScalar(*v); v++; 1371 } 1372 for (j=0; j<A->n; j++) { 1373 if (tmp[j] > *nrm) *nrm = tmp[j]; 1374 } 1375 ierr = PetscFree(tmp);CHKERRQ(ierr); 1376 } else if (type == NORM_INFINITY) { 1377 *nrm = 0.0; 1378 for (j=0; j<A->m; j++) { 1379 v = a->a + a->i[j]; 1380 sum = 0.0; 1381 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1382 sum += PetscAbsScalar(*v); v++; 1383 } 1384 if (sum > *nrm) *nrm = sum; 1385 } 1386 } else { 1387 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1388 } 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "MatTranspose_SeqAIJ" 1394 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1395 { 1396 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1397 Mat C; 1398 PetscErrorCode ierr; 1399 PetscInt i,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1400 PetscScalar *array = a->a; 1401 1402 PetscFunctionBegin; 1403 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1404 ierr = PetscMalloc((1+A->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1405 ierr = PetscMemzero(col,(1+A->n)*sizeof(PetscInt));CHKERRQ(ierr); 1406 1407 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1408 ierr = MatCreate(A->comm,A->n,m,A->n,m,&C);CHKERRQ(ierr); 1409 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1410 ierr = MatSeqAIJSetPreallocation(C,0,col);CHKERRQ(ierr); 1411 ierr = PetscFree(col);CHKERRQ(ierr); 1412 for (i=0; i<m; i++) { 1413 len = ai[i+1]-ai[i]; 1414 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1415 array += len; 1416 aj += len; 1417 } 1418 1419 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1420 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1421 1422 if (B) { 1423 *B = C; 1424 } else { 1425 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1426 } 1427 PetscFunctionReturn(0); 1428 } 1429 1430 EXTERN_C_BEGIN 1431 #undef __FUNCT__ 1432 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1433 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1434 { 1435 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1436 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1437 PetscErrorCode ierr; 1438 PetscInt ma,na,mb,nb, i; 1439 1440 PetscFunctionBegin; 1441 bij = (Mat_SeqAIJ *) B->data; 1442 1443 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1444 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1445 if (ma!=nb || na!=mb){ 1446 *f = PETSC_FALSE; 1447 PetscFunctionReturn(0); 1448 } 1449 aii = aij->i; bii = bij->i; 1450 adx = aij->j; bdx = bij->j; 1451 va = aij->a; vb = bij->a; 1452 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1453 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1454 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1455 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1456 1457 *f = PETSC_TRUE; 1458 for (i=0; i<ma; i++) { 1459 while (aptr[i]<aii[i+1]) { 1460 PetscInt idc,idr; 1461 PetscScalar vc,vr; 1462 /* column/row index/value */ 1463 idc = adx[aptr[i]]; 1464 idr = bdx[bptr[idc]]; 1465 vc = va[aptr[i]]; 1466 vr = vb[bptr[idc]]; 1467 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1468 *f = PETSC_FALSE; 1469 goto done; 1470 } else { 1471 aptr[i]++; 1472 if (B || i!=idc) bptr[idc]++; 1473 } 1474 } 1475 } 1476 done: 1477 ierr = PetscFree(aptr);CHKERRQ(ierr); 1478 if (B) { 1479 ierr = PetscFree(bptr);CHKERRQ(ierr); 1480 } 1481 PetscFunctionReturn(0); 1482 } 1483 EXTERN_C_END 1484 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1487 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1488 { 1489 PetscErrorCode ierr; 1490 PetscFunctionBegin; 1491 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1497 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1498 { 1499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1500 PetscScalar *l,*r,x,*v; 1501 PetscErrorCode ierr; 1502 PetscInt i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1503 1504 PetscFunctionBegin; 1505 if (ll) { 1506 /* The local size is used so that VecMPI can be passed to this routine 1507 by MatDiagonalScale_MPIAIJ */ 1508 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1509 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1510 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1511 v = a->a; 1512 for (i=0; i<m; i++) { 1513 x = l[i]; 1514 M = a->i[i+1] - a->i[i]; 1515 for (j=0; j<M; j++) { (*v++) *= x;} 1516 } 1517 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1518 PetscLogFlops(nz); 1519 } 1520 if (rr) { 1521 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1522 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1523 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1524 v = a->a; jj = a->j; 1525 for (i=0; i<nz; i++) { 1526 (*v++) *= r[*jj++]; 1527 } 1528 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1529 PetscLogFlops(nz); 1530 } 1531 PetscFunctionReturn(0); 1532 } 1533 1534 #undef __FUNCT__ 1535 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1536 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1537 { 1538 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1539 PetscErrorCode ierr; 1540 PetscInt *smap,i,k,kstart,kend,oldcols = A->n,*lens; 1541 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1542 PetscInt *irow,*icol,nrows,ncols; 1543 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1544 PetscScalar *a_new,*mat_a; 1545 Mat C; 1546 PetscTruth stride; 1547 1548 PetscFunctionBegin; 1549 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1550 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1551 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1552 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1553 1554 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1555 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1556 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1557 1558 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1559 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1560 if (stride && step == 1) { 1561 /* special case of contiguous rows */ 1562 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1563 starts = lens + nrows; 1564 /* loop over new rows determining lens and starting points */ 1565 for (i=0; i<nrows; i++) { 1566 kstart = ai[irow[i]]; 1567 kend = kstart + ailen[irow[i]]; 1568 for (k=kstart; k<kend; k++) { 1569 if (aj[k] >= first) { 1570 starts[i] = k; 1571 break; 1572 } 1573 } 1574 sum = 0; 1575 while (k < kend) { 1576 if (aj[k++] >= first+ncols) break; 1577 sum++; 1578 } 1579 lens[i] = sum; 1580 } 1581 /* create submatrix */ 1582 if (scall == MAT_REUSE_MATRIX) { 1583 PetscInt n_cols,n_rows; 1584 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1585 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1586 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1587 C = *B; 1588 } else { 1589 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1590 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1591 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1592 } 1593 c = (Mat_SeqAIJ*)C->data; 1594 1595 /* loop over rows inserting into submatrix */ 1596 a_new = c->a; 1597 j_new = c->j; 1598 i_new = c->i; 1599 1600 for (i=0; i<nrows; i++) { 1601 ii = starts[i]; 1602 lensi = lens[i]; 1603 for (k=0; k<lensi; k++) { 1604 *j_new++ = aj[ii+k] - first; 1605 } 1606 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1607 a_new += lensi; 1608 i_new[i+1] = i_new[i] + lensi; 1609 c->ilen[i] = lensi; 1610 } 1611 ierr = PetscFree(lens);CHKERRQ(ierr); 1612 } else { 1613 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1614 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1615 1616 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1617 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1618 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1619 /* determine lens of each row */ 1620 for (i=0; i<nrows; i++) { 1621 kstart = ai[irow[i]]; 1622 kend = kstart + a->ilen[irow[i]]; 1623 lens[i] = 0; 1624 for (k=kstart; k<kend; k++) { 1625 if (smap[aj[k]]) { 1626 lens[i]++; 1627 } 1628 } 1629 } 1630 /* Create and fill new matrix */ 1631 if (scall == MAT_REUSE_MATRIX) { 1632 PetscTruth equal; 1633 1634 c = (Mat_SeqAIJ *)((*B)->data); 1635 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1636 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1637 if (!equal) { 1638 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1639 } 1640 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(PetscInt));CHKERRQ(ierr); 1641 C = *B; 1642 } else { 1643 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1644 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1645 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1646 } 1647 c = (Mat_SeqAIJ *)(C->data); 1648 for (i=0; i<nrows; i++) { 1649 row = irow[i]; 1650 kstart = ai[row]; 1651 kend = kstart + a->ilen[row]; 1652 mat_i = c->i[i]; 1653 mat_j = c->j + mat_i; 1654 mat_a = c->a + mat_i; 1655 mat_ilen = c->ilen + i; 1656 for (k=kstart; k<kend; k++) { 1657 if ((tcol=smap[a->j[k]])) { 1658 *mat_j++ = tcol - 1; 1659 *mat_a++ = a->a[k]; 1660 (*mat_ilen)++; 1661 1662 } 1663 } 1664 } 1665 /* Free work space */ 1666 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1667 ierr = PetscFree(smap);CHKERRQ(ierr); 1668 ierr = PetscFree(lens);CHKERRQ(ierr); 1669 } 1670 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1671 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1672 1673 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1674 *B = C; 1675 PetscFunctionReturn(0); 1676 } 1677 1678 /* 1679 */ 1680 #undef __FUNCT__ 1681 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1682 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1683 { 1684 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1685 PetscErrorCode ierr; 1686 Mat outA; 1687 PetscTruth row_identity,col_identity; 1688 1689 PetscFunctionBegin; 1690 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1691 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1692 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1693 if (!row_identity || !col_identity) { 1694 SETERRQ(PETSC_ERR_ARG_WRONG,"Row and column permutations must be identity for in-place ILU"); 1695 } 1696 1697 outA = inA; 1698 inA->factor = FACTOR_LU; 1699 a->row = row; 1700 a->col = col; 1701 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1702 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1703 1704 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1705 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1706 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1707 PetscLogObjectParent(inA,a->icol); 1708 1709 if (!a->solve_work) { /* this matrix may have been factored before */ 1710 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1711 } 1712 1713 if (!a->diag) { 1714 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1715 } 1716 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1717 PetscFunctionReturn(0); 1718 } 1719 1720 #include "petscblaslapack.h" 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatScale_SeqAIJ" 1723 PetscErrorCode MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1724 { 1725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1726 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1727 1728 PetscFunctionBegin; 1729 BLscal_(&bnz,(PetscScalar*)alpha,a->a,&one); 1730 PetscLogFlops(a->nz); 1731 PetscFunctionReturn(0); 1732 } 1733 1734 #undef __FUNCT__ 1735 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1736 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1737 { 1738 PetscErrorCode ierr; 1739 PetscInt i; 1740 1741 PetscFunctionBegin; 1742 if (scall == MAT_INITIAL_MATRIX) { 1743 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1744 } 1745 1746 for (i=0; i<n; i++) { 1747 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1748 } 1749 PetscFunctionReturn(0); 1750 } 1751 1752 #undef __FUNCT__ 1753 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1754 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1755 { 1756 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1757 PetscErrorCode ierr; 1758 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1759 PetscInt start,end,*ai,*aj; 1760 PetscBT table; 1761 1762 PetscFunctionBegin; 1763 m = A->m; 1764 ai = a->i; 1765 aj = a->j; 1766 1767 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1768 1769 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1770 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1771 1772 for (i=0; i<is_max; i++) { 1773 /* Initialize the two local arrays */ 1774 isz = 0; 1775 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1776 1777 /* Extract the indices, assume there can be duplicate entries */ 1778 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1779 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1780 1781 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1782 for (j=0; j<n ; ++j){ 1783 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1784 } 1785 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1786 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1787 1788 k = 0; 1789 for (j=0; j<ov; j++){ /* for each overlap */ 1790 n = isz; 1791 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1792 row = nidx[k]; 1793 start = ai[row]; 1794 end = ai[row+1]; 1795 for (l = start; l<end ; l++){ 1796 val = aj[l] ; 1797 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1798 } 1799 } 1800 } 1801 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1802 } 1803 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1804 ierr = PetscFree(nidx);CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 /* -------------------------------------------------------------- */ 1809 #undef __FUNCT__ 1810 #define __FUNCT__ "MatPermute_SeqAIJ" 1811 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1812 { 1813 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1814 PetscErrorCode ierr; 1815 PetscInt i,nz,m = A->m,n = A->n,*col; 1816 PetscInt *row,*cnew,j,*lens; 1817 IS icolp,irowp; 1818 PetscInt *cwork; 1819 PetscScalar *vwork; 1820 1821 PetscFunctionBegin; 1822 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1823 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1824 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1825 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1826 1827 /* determine lengths of permuted rows */ 1828 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1829 for (i=0; i<m; i++) { 1830 lens[row[i]] = a->i[i+1] - a->i[i]; 1831 } 1832 ierr = MatCreate(A->comm,m,n,m,n,B);CHKERRQ(ierr); 1833 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1834 ierr = MatSeqAIJSetPreallocation(*B,0,lens);CHKERRQ(ierr); 1835 ierr = PetscFree(lens);CHKERRQ(ierr); 1836 1837 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1838 for (i=0; i<m; i++) { 1839 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1840 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1841 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1842 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1843 } 1844 ierr = PetscFree(cnew);CHKERRQ(ierr); 1845 (*B)->assembled = PETSC_FALSE; 1846 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1847 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1848 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1849 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1850 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1851 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1857 PetscErrorCode MatPrintHelp_SeqAIJ(Mat A) 1858 { 1859 static PetscTruth called = PETSC_FALSE; 1860 MPI_Comm comm = A->comm; 1861 PetscErrorCode ierr; 1862 1863 PetscFunctionBegin; 1864 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1865 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1866 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1867 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1868 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1869 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1870 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_compressedrow: Do not use compressedrow\n");CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatCopy_SeqAIJ" 1876 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1877 { 1878 PetscErrorCode ierr; 1879 1880 PetscFunctionBegin; 1881 /* If the two matrices have the same copy implementation, use fast copy. */ 1882 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1883 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1884 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1885 1886 if (a->i[A->m] != b->i[B->m]) { 1887 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1888 } 1889 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1890 } else { 1891 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1892 } 1893 PetscFunctionReturn(0); 1894 } 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1898 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1899 { 1900 PetscErrorCode ierr; 1901 1902 PetscFunctionBegin; 1903 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1904 PetscFunctionReturn(0); 1905 } 1906 1907 #undef __FUNCT__ 1908 #define __FUNCT__ "MatGetArray_SeqAIJ" 1909 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1910 { 1911 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1912 PetscFunctionBegin; 1913 *array = a->a; 1914 PetscFunctionReturn(0); 1915 } 1916 1917 #undef __FUNCT__ 1918 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1919 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1920 { 1921 PetscFunctionBegin; 1922 PetscFunctionReturn(0); 1923 } 1924 1925 #undef __FUNCT__ 1926 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1927 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1928 { 1929 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 1930 PetscErrorCode ierr; 1931 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1932 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1933 PetscScalar *vscale_array; 1934 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1935 Vec w1,w2,w3; 1936 void *fctx = coloring->fctx; 1937 PetscTruth flg; 1938 1939 PetscFunctionBegin; 1940 if (!coloring->w1) { 1941 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1942 PetscLogObjectParent(coloring,coloring->w1); 1943 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1944 PetscLogObjectParent(coloring,coloring->w2); 1945 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1946 PetscLogObjectParent(coloring,coloring->w3); 1947 } 1948 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1949 1950 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1951 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1952 if (flg) { 1953 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1954 } else { 1955 PetscTruth assembled; 1956 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 1957 if (assembled) { 1958 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1959 } 1960 } 1961 1962 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1963 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1964 1965 /* 1966 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1967 coloring->F for the coarser grids from the finest 1968 */ 1969 if (coloring->F) { 1970 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1971 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1972 if (m1 != m2) { 1973 coloring->F = 0; 1974 } 1975 } 1976 1977 if (coloring->F) { 1978 w1 = coloring->F; 1979 coloring->F = 0; 1980 } else { 1981 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1982 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1983 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1984 } 1985 1986 /* 1987 Compute all the scale factors and share with other processors 1988 */ 1989 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1990 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1991 for (k=0; k<coloring->ncolors; k++) { 1992 /* 1993 Loop over each column associated with color adding the 1994 perturbation to the vector w3. 1995 */ 1996 for (l=0; l<coloring->ncolumns[k]; l++) { 1997 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1998 dx = xx[col]; 1999 if (dx == 0.0) dx = 1.0; 2000 #if !defined(PETSC_USE_COMPLEX) 2001 if (dx < umin && dx >= 0.0) dx = umin; 2002 else if (dx < 0.0 && dx > -umin) dx = -umin; 2003 #else 2004 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2005 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2006 #endif 2007 dx *= epsilon; 2008 vscale_array[col] = 1.0/dx; 2009 } 2010 } 2011 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2012 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2013 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2014 2015 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2016 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2017 2018 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2019 else vscaleforrow = coloring->columnsforrow; 2020 2021 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2022 /* 2023 Loop over each color 2024 */ 2025 for (k=0; k<coloring->ncolors; k++) { 2026 coloring->currentcolor = k; 2027 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2028 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2029 /* 2030 Loop over each column associated with color adding the 2031 perturbation to the vector w3. 2032 */ 2033 for (l=0; l<coloring->ncolumns[k]; l++) { 2034 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2035 dx = xx[col]; 2036 if (dx == 0.0) dx = 1.0; 2037 #if !defined(PETSC_USE_COMPLEX) 2038 if (dx < umin && dx >= 0.0) dx = umin; 2039 else if (dx < 0.0 && dx > -umin) dx = -umin; 2040 #else 2041 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2042 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2043 #endif 2044 dx *= epsilon; 2045 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2046 w3_array[col] += dx; 2047 } 2048 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2049 2050 /* 2051 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2052 */ 2053 2054 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2055 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2056 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2057 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2058 2059 /* 2060 Loop over rows of vector, putting results into Jacobian matrix 2061 */ 2062 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2063 for (l=0; l<coloring->nrows[k]; l++) { 2064 row = coloring->rows[k][l]; 2065 col = coloring->columnsforrow[k][l]; 2066 y[row] *= vscale_array[vscaleforrow[k][l]]; 2067 srow = row + start; 2068 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2069 } 2070 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2071 } 2072 coloring->currentcolor = k; 2073 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2074 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2075 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2076 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2077 PetscFunctionReturn(0); 2078 } 2079 2080 #include "petscblaslapack.h" 2081 #undef __FUNCT__ 2082 #define __FUNCT__ "MatAXPY_SeqAIJ" 2083 PetscErrorCode MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2084 { 2085 PetscErrorCode ierr; 2086 PetscInt i; 2087 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2088 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2089 2090 PetscFunctionBegin; 2091 if (str == SAME_NONZERO_PATTERN) { 2092 BLaxpy_(&bnz,(PetscScalar*)a,x->a,&one,y->a,&one); 2093 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2094 if (y->xtoy && y->XtoY != X) { 2095 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2096 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2097 } 2098 if (!y->xtoy) { /* get xtoy */ 2099 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2100 y->XtoY = X; 2101 } 2102 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2103 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2104 } else { 2105 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2106 } 2107 PetscFunctionReturn(0); 2108 } 2109 2110 #undef __FUNCT__ 2111 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2112 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2113 { 2114 PetscFunctionBegin; 2115 PetscFunctionReturn(0); 2116 } 2117 2118 /* -------------------------------------------------------------------*/ 2119 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2120 MatGetRow_SeqAIJ, 2121 MatRestoreRow_SeqAIJ, 2122 MatMult_SeqAIJ, 2123 /* 4*/ MatMultAdd_SeqAIJ, 2124 MatMultTranspose_SeqAIJ, 2125 MatMultTransposeAdd_SeqAIJ, 2126 MatSolve_SeqAIJ, 2127 MatSolveAdd_SeqAIJ, 2128 MatSolveTranspose_SeqAIJ, 2129 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2130 MatLUFactor_SeqAIJ, 2131 0, 2132 MatRelax_SeqAIJ, 2133 MatTranspose_SeqAIJ, 2134 /*15*/ MatGetInfo_SeqAIJ, 2135 MatEqual_SeqAIJ, 2136 MatGetDiagonal_SeqAIJ, 2137 MatDiagonalScale_SeqAIJ, 2138 MatNorm_SeqAIJ, 2139 /*20*/ 0, 2140 MatAssemblyEnd_SeqAIJ, 2141 MatCompress_SeqAIJ, 2142 MatSetOption_SeqAIJ, 2143 MatZeroEntries_SeqAIJ, 2144 /*25*/ MatZeroRows_SeqAIJ, 2145 MatLUFactorSymbolic_SeqAIJ, 2146 MatLUFactorNumeric_SeqAIJ, 2147 MatCholeskyFactorSymbolic_SeqAIJ, 2148 MatCholeskyFactorNumeric_SeqAIJ, 2149 /*30*/ MatSetUpPreallocation_SeqAIJ, 2150 MatILUFactorSymbolic_SeqAIJ, 2151 MatICCFactorSymbolic_SeqAIJ, 2152 MatGetArray_SeqAIJ, 2153 MatRestoreArray_SeqAIJ, 2154 /*35*/ MatDuplicate_SeqAIJ, 2155 0, 2156 0, 2157 MatILUFactor_SeqAIJ, 2158 0, 2159 /*40*/ MatAXPY_SeqAIJ, 2160 MatGetSubMatrices_SeqAIJ, 2161 MatIncreaseOverlap_SeqAIJ, 2162 MatGetValues_SeqAIJ, 2163 MatCopy_SeqAIJ, 2164 /*45*/ MatPrintHelp_SeqAIJ, 2165 MatScale_SeqAIJ, 2166 0, 2167 0, 2168 MatILUDTFactor_SeqAIJ, 2169 /*50*/ MatSetBlockSize_SeqAIJ, 2170 MatGetRowIJ_SeqAIJ, 2171 MatRestoreRowIJ_SeqAIJ, 2172 MatGetColumnIJ_SeqAIJ, 2173 MatRestoreColumnIJ_SeqAIJ, 2174 /*55*/ MatFDColoringCreate_SeqAIJ, 2175 0, 2176 0, 2177 MatPermute_SeqAIJ, 2178 0, 2179 /*60*/ 0, 2180 MatDestroy_SeqAIJ, 2181 MatView_SeqAIJ, 2182 MatGetPetscMaps_Petsc, 2183 0, 2184 /*65*/ 0, 2185 0, 2186 0, 2187 0, 2188 0, 2189 /*70*/ 0, 2190 0, 2191 MatSetColoring_SeqAIJ, 2192 MatSetValuesAdic_SeqAIJ, 2193 MatSetValuesAdifor_SeqAIJ, 2194 /*75*/ MatFDColoringApply_SeqAIJ, 2195 0, 2196 0, 2197 0, 2198 0, 2199 /*80*/ 0, 2200 0, 2201 0, 2202 0, 2203 MatLoad_SeqAIJ, 2204 /*85*/ MatIsSymmetric_SeqAIJ, 2205 0, 2206 0, 2207 0, 2208 0, 2209 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2210 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2211 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2212 MatPtAP_SeqAIJ_SeqAIJ, 2213 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2214 /*95*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2215 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2216 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2217 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2218 }; 2219 2220 EXTERN_C_BEGIN 2221 #undef __FUNCT__ 2222 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2223 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2224 { 2225 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2226 PetscInt i,nz,n; 2227 2228 PetscFunctionBegin; 2229 2230 nz = aij->maxnz; 2231 n = mat->n; 2232 for (i=0; i<nz; i++) { 2233 aij->j[i] = indices[i]; 2234 } 2235 aij->nz = nz; 2236 for (i=0; i<n; i++) { 2237 aij->ilen[i] = aij->imax[i]; 2238 } 2239 2240 PetscFunctionReturn(0); 2241 } 2242 EXTERN_C_END 2243 2244 #undef __FUNCT__ 2245 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2246 /*@ 2247 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2248 in the matrix. 2249 2250 Input Parameters: 2251 + mat - the SeqAIJ matrix 2252 - indices - the column indices 2253 2254 Level: advanced 2255 2256 Notes: 2257 This can be called if you have precomputed the nonzero structure of the 2258 matrix and want to provide it to the matrix object to improve the performance 2259 of the MatSetValues() operation. 2260 2261 You MUST have set the correct numbers of nonzeros per row in the call to 2262 MatCreateSeqAIJ(). 2263 2264 MUST be called before any calls to MatSetValues(); 2265 2266 The indices should start with zero, not one. 2267 2268 @*/ 2269 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2270 { 2271 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2272 2273 PetscFunctionBegin; 2274 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2275 PetscValidPointer(indices,2); 2276 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2277 if (f) { 2278 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2279 } else { 2280 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2281 } 2282 PetscFunctionReturn(0); 2283 } 2284 2285 /* ----------------------------------------------------------------------------------------*/ 2286 2287 EXTERN_C_BEGIN 2288 #undef __FUNCT__ 2289 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2290 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 2291 { 2292 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2293 PetscErrorCode ierr; 2294 size_t nz = aij->i[mat->m]; 2295 2296 PetscFunctionBegin; 2297 if (aij->nonew != 1) { 2298 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2299 } 2300 2301 /* allocate space for values if not already there */ 2302 if (!aij->saved_values) { 2303 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2304 } 2305 2306 /* copy values over */ 2307 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2308 PetscFunctionReturn(0); 2309 } 2310 EXTERN_C_END 2311 2312 #undef __FUNCT__ 2313 #define __FUNCT__ "MatStoreValues" 2314 /*@ 2315 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2316 example, reuse of the linear part of a Jacobian, while recomputing the 2317 nonlinear portion. 2318 2319 Collect on Mat 2320 2321 Input Parameters: 2322 . mat - the matrix (currently on AIJ matrices support this option) 2323 2324 Level: advanced 2325 2326 Common Usage, with SNESSolve(): 2327 $ Create Jacobian matrix 2328 $ Set linear terms into matrix 2329 $ Apply boundary conditions to matrix, at this time matrix must have 2330 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2331 $ boundary conditions again will not change the nonzero structure 2332 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2333 $ ierr = MatStoreValues(mat); 2334 $ Call SNESSetJacobian() with matrix 2335 $ In your Jacobian routine 2336 $ ierr = MatRetrieveValues(mat); 2337 $ Set nonlinear terms in matrix 2338 2339 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2340 $ // build linear portion of Jacobian 2341 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2342 $ ierr = MatStoreValues(mat); 2343 $ loop over nonlinear iterations 2344 $ ierr = MatRetrieveValues(mat); 2345 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2346 $ // call MatAssemblyBegin/End() on matrix 2347 $ Solve linear system with Jacobian 2348 $ endloop 2349 2350 Notes: 2351 Matrix must already be assemblied before calling this routine 2352 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2353 calling this routine. 2354 2355 .seealso: MatRetrieveValues() 2356 2357 @*/ 2358 PetscErrorCode MatStoreValues(Mat mat) 2359 { 2360 PetscErrorCode ierr,(*f)(Mat); 2361 2362 PetscFunctionBegin; 2363 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2364 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2365 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2366 2367 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2368 if (f) { 2369 ierr = (*f)(mat);CHKERRQ(ierr); 2370 } else { 2371 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2372 } 2373 PetscFunctionReturn(0); 2374 } 2375 2376 EXTERN_C_BEGIN 2377 #undef __FUNCT__ 2378 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2379 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 2380 { 2381 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2382 PetscErrorCode ierr; 2383 PetscInt nz = aij->i[mat->m]; 2384 2385 PetscFunctionBegin; 2386 if (aij->nonew != 1) { 2387 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2388 } 2389 if (!aij->saved_values) { 2390 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2391 } 2392 /* copy values over */ 2393 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396 EXTERN_C_END 2397 2398 #undef __FUNCT__ 2399 #define __FUNCT__ "MatRetrieveValues" 2400 /*@ 2401 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2402 example, reuse of the linear part of a Jacobian, while recomputing the 2403 nonlinear portion. 2404 2405 Collect on Mat 2406 2407 Input Parameters: 2408 . mat - the matrix (currently on AIJ matrices support this option) 2409 2410 Level: advanced 2411 2412 .seealso: MatStoreValues() 2413 2414 @*/ 2415 PetscErrorCode MatRetrieveValues(Mat mat) 2416 { 2417 PetscErrorCode ierr,(*f)(Mat); 2418 2419 PetscFunctionBegin; 2420 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2421 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2422 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2423 2424 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2425 if (f) { 2426 ierr = (*f)(mat);CHKERRQ(ierr); 2427 } else { 2428 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2429 } 2430 PetscFunctionReturn(0); 2431 } 2432 2433 2434 /* --------------------------------------------------------------------------------*/ 2435 #undef __FUNCT__ 2436 #define __FUNCT__ "MatCreateSeqAIJ" 2437 /*@C 2438 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2439 (the default parallel PETSc format). For good matrix assembly performance 2440 the user should preallocate the matrix storage by setting the parameter nz 2441 (or the array nnz). By setting these parameters accurately, performance 2442 during matrix assembly can be increased by more than a factor of 50. 2443 2444 Collective on MPI_Comm 2445 2446 Input Parameters: 2447 + comm - MPI communicator, set to PETSC_COMM_SELF 2448 . m - number of rows 2449 . n - number of columns 2450 . nz - number of nonzeros per row (same for all rows) 2451 - nnz - array containing the number of nonzeros in the various rows 2452 (possibly different for each row) or PETSC_NULL 2453 2454 Output Parameter: 2455 . A - the matrix 2456 2457 Notes: 2458 If nnz is given then nz is ignored 2459 2460 The AIJ format (also called the Yale sparse matrix format or 2461 compressed row storage), is fully compatible with standard Fortran 77 2462 storage. That is, the stored row and column indices can begin at 2463 either one (as in Fortran) or zero. See the users' manual for details. 2464 2465 Specify the preallocated storage with either nz or nnz (not both). 2466 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2467 allocation. For large problems you MUST preallocate memory or you 2468 will get TERRIBLE performance, see the users' manual chapter on matrices. 2469 2470 By default, this format uses inodes (identical nodes) when possible, to 2471 improve numerical efficiency of matrix-vector products and solves. We 2472 search for consecutive rows with the same nonzero structure, thereby 2473 reusing matrix information to achieve increased efficiency. 2474 2475 Options Database Keys: 2476 + -mat_aij_no_inode - Do not use inodes 2477 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2478 - -mat_aij_oneindex - Internally use indexing starting at 1 2479 rather than 0. Note that when calling MatSetValues(), 2480 the user still MUST index entries starting at 0! 2481 2482 Level: intermediate 2483 2484 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2485 2486 @*/ 2487 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2488 { 2489 PetscErrorCode ierr; 2490 2491 PetscFunctionBegin; 2492 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2493 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2494 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 #define SKIP_ALLOCATION -4 2499 2500 #undef __FUNCT__ 2501 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2502 /*@C 2503 MatSeqAIJSetPreallocation - For good matrix assembly performance 2504 the user should preallocate the matrix storage by setting the parameter nz 2505 (or the array nnz). By setting these parameters accurately, performance 2506 during matrix assembly can be increased by more than a factor of 50. 2507 2508 Collective on MPI_Comm 2509 2510 Input Parameters: 2511 + comm - MPI communicator, set to PETSC_COMM_SELF 2512 . m - number of rows 2513 . n - number of columns 2514 . nz - number of nonzeros per row (same for all rows) 2515 - nnz - array containing the number of nonzeros in the various rows 2516 (possibly different for each row) or PETSC_NULL 2517 2518 Output Parameter: 2519 . A - the matrix 2520 2521 Notes: 2522 If nnz is given then nz is ignored 2523 2524 The AIJ format (also called the Yale sparse matrix format or 2525 compressed row storage), is fully compatible with standard Fortran 77 2526 storage. That is, the stored row and column indices can begin at 2527 either one (as in Fortran) or zero. See the users' manual for details. 2528 2529 Specify the preallocated storage with either nz or nnz (not both). 2530 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2531 allocation. For large problems you MUST preallocate memory or you 2532 will get TERRIBLE performance, see the users' manual chapter on matrices. 2533 2534 By default, this format uses inodes (identical nodes) when possible, to 2535 improve numerical efficiency of matrix-vector products and solves. We 2536 search for consecutive rows with the same nonzero structure, thereby 2537 reusing matrix information to achieve increased efficiency. 2538 2539 Options Database Keys: 2540 + -mat_aij_no_inode - Do not use inodes 2541 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2542 - -mat_aij_oneindex - Internally use indexing starting at 1 2543 rather than 0. Note that when calling MatSetValues(), 2544 the user still MUST index entries starting at 0! 2545 2546 Level: intermediate 2547 2548 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2549 2550 @*/ 2551 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2552 { 2553 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2554 2555 PetscFunctionBegin; 2556 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2557 if (f) { 2558 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2559 } 2560 PetscFunctionReturn(0); 2561 } 2562 2563 EXTERN_C_BEGIN 2564 #undef __FUNCT__ 2565 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2566 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2567 { 2568 Mat_SeqAIJ *b; 2569 size_t len = 0; 2570 PetscTruth skipallocation = PETSC_FALSE; 2571 PetscErrorCode ierr; 2572 PetscInt i; 2573 2574 PetscFunctionBegin; 2575 2576 if (nz == SKIP_ALLOCATION) { 2577 skipallocation = PETSC_TRUE; 2578 nz = 0; 2579 } 2580 2581 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2582 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2583 if (nnz) { 2584 for (i=0; i<B->m; i++) { 2585 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2586 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2587 } 2588 } 2589 2590 B->preallocated = PETSC_TRUE; 2591 b = (Mat_SeqAIJ*)B->data; 2592 2593 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->imax);CHKERRQ(ierr); 2594 if (!nnz) { 2595 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2596 else if (nz <= 0) nz = 1; 2597 for (i=0; i<B->m; i++) b->imax[i] = nz; 2598 nz = nz*B->m; 2599 } else { 2600 nz = 0; 2601 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2602 } 2603 2604 if (!skipallocation) { 2605 /* allocate the matrix space */ 2606 len = ((size_t) nz)*(sizeof(PetscInt) + sizeof(PetscScalar)) + (B->m+1)*sizeof(PetscInt); 2607 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2608 b->j = (PetscInt*)(b->a + nz); 2609 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2610 b->i = b->j + nz; 2611 b->i[0] = 0; 2612 for (i=1; i<B->m+1; i++) { 2613 b->i[i] = b->i[i-1] + b->imax[i-1]; 2614 } 2615 b->singlemalloc = PETSC_TRUE; 2616 b->freedata = PETSC_TRUE; 2617 } else { 2618 b->freedata = PETSC_FALSE; 2619 } 2620 2621 /* b->ilen will count nonzeros in each row so far. */ 2622 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->ilen);CHKERRQ(ierr); 2623 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2624 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2625 2626 b->nz = 0; 2627 b->maxnz = nz; 2628 B->info.nz_unneeded = (double)b->maxnz; 2629 PetscFunctionReturn(0); 2630 } 2631 EXTERN_C_END 2632 2633 /*MC 2634 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2635 based on compressed sparse row format. 2636 2637 Options Database Keys: 2638 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2639 2640 Level: beginner 2641 2642 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2643 M*/ 2644 2645 EXTERN_C_BEGIN 2646 #undef __FUNCT__ 2647 #define __FUNCT__ "MatCreate_SeqAIJ" 2648 PetscErrorCode MatCreate_SeqAIJ(Mat B) 2649 { 2650 Mat_SeqAIJ *b; 2651 PetscErrorCode ierr; 2652 PetscMPIInt size; 2653 2654 PetscFunctionBegin; 2655 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2656 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2657 2658 B->m = B->M = PetscMax(B->m,B->M); 2659 B->n = B->N = PetscMax(B->n,B->N); 2660 2661 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2662 B->data = (void*)b; 2663 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2664 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2665 B->factor = 0; 2666 B->lupivotthreshold = 1.0; 2667 B->mapping = 0; 2668 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2669 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2670 b->row = 0; 2671 b->col = 0; 2672 b->icol = 0; 2673 b->reallocs = 0; 2674 b->lu_shift = PETSC_FALSE; 2675 2676 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2677 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2678 2679 b->sorted = PETSC_FALSE; 2680 b->ignorezeroentries = PETSC_FALSE; 2681 b->roworiented = PETSC_TRUE; 2682 b->nonew = 0; 2683 b->diag = 0; 2684 b->solve_work = 0; 2685 B->spptr = 0; 2686 b->inode.use = PETSC_TRUE; 2687 b->inode.node_count = 0; 2688 b->inode.size = 0; 2689 b->inode.limit = 5; 2690 b->inode.max_limit = 5; 2691 b->saved_values = 0; 2692 b->idiag = 0; 2693 b->ssor = 0; 2694 b->keepzeroedrows = PETSC_FALSE; 2695 b->xtoy = 0; 2696 b->XtoY = 0; 2697 b->compressedrow.use = PETSC_TRUE; 2698 b->compressedrow.nrows = B->m; 2699 b->compressedrow.i = PETSC_NULL; 2700 b->compressedrow.rindex = PETSC_NULL; 2701 b->compressedrow.checked = PETSC_FALSE; 2702 2703 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2704 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2706 "MatSeqAIJSetColumnIndices_SeqAIJ", 2707 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2709 "MatStoreValues_SeqAIJ", 2710 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2711 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2712 "MatRetrieveValues_SeqAIJ", 2713 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2714 #if !defined(PETSC_USE_64BIT_INT) 2715 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2716 "MatConvert_SeqAIJ_SeqSBAIJ", 2717 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2718 #endif 2719 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2720 "MatConvert_SeqAIJ_SeqBAIJ", 2721 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2722 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2723 "MatIsTranspose_SeqAIJ", 2724 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2725 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2726 "MatSeqAIJSetPreallocation_SeqAIJ", 2727 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2728 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2729 "MatReorderForNonzeroDiagonal_SeqAIJ", 2730 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2732 "MatAdjustForInodes_SeqAIJ", 2733 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2735 "MatSeqAIJGetInodeSizes_SeqAIJ", 2736 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2737 PetscFunctionReturn(0); 2738 } 2739 EXTERN_C_END 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2743 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2744 { 2745 Mat C; 2746 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2747 PetscErrorCode ierr; 2748 PetscInt i,m = A->m; 2749 size_t len; 2750 2751 PetscFunctionBegin; 2752 *B = 0; 2753 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2754 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2755 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2756 2757 c = (Mat_SeqAIJ*)C->data; 2758 2759 C->factor = A->factor; 2760 c->row = 0; 2761 c->col = 0; 2762 c->icol = 0; 2763 c->keepzeroedrows = a->keepzeroedrows; 2764 C->assembled = PETSC_TRUE; 2765 2766 C->M = A->m; 2767 C->N = A->n; 2768 C->bs = A->bs; 2769 2770 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->imax);CHKERRQ(ierr); 2771 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->ilen);CHKERRQ(ierr); 2772 for (i=0; i<m; i++) { 2773 c->imax[i] = a->imax[i]; 2774 c->ilen[i] = a->ilen[i]; 2775 } 2776 2777 /* allocate the matrix space */ 2778 c->singlemalloc = PETSC_TRUE; 2779 len = ((size_t) (m+1))*sizeof(PetscInt)+(a->i[m])*(sizeof(PetscScalar)+sizeof(PetscInt)); 2780 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2781 c->j = (PetscInt*)(c->a + a->i[m] ); 2782 c->i = c->j + a->i[m]; 2783 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2784 if (m > 0) { 2785 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2786 if (cpvalues == MAT_COPY_VALUES) { 2787 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2788 } else { 2789 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2790 } 2791 } 2792 2793 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2794 c->sorted = a->sorted; 2795 c->roworiented = a->roworiented; 2796 c->nonew = a->nonew; 2797 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2798 c->saved_values = 0; 2799 c->idiag = 0; 2800 c->ssor = 0; 2801 c->ignorezeroentries = a->ignorezeroentries; 2802 c->freedata = PETSC_TRUE; 2803 2804 if (a->diag) { 2805 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2806 PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt)); 2807 for (i=0; i<m; i++) { 2808 c->diag[i] = a->diag[i]; 2809 } 2810 } else c->diag = 0; 2811 c->inode.use = a->inode.use; 2812 c->inode.limit = a->inode.limit; 2813 c->inode.max_limit = a->inode.max_limit; 2814 if (a->inode.size){ 2815 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->inode.size);CHKERRQ(ierr); 2816 c->inode.node_count = a->inode.node_count; 2817 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2818 } else { 2819 c->inode.size = 0; 2820 c->inode.node_count = 0; 2821 } 2822 c->nz = a->nz; 2823 c->maxnz = a->maxnz; 2824 c->solve_work = 0; 2825 C->preallocated = PETSC_TRUE; 2826 2827 *B = C; 2828 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2829 PetscFunctionReturn(0); 2830 } 2831 2832 #undef __FUNCT__ 2833 #define __FUNCT__ "MatLoad_SeqAIJ" 2834 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2835 { 2836 Mat_SeqAIJ *a; 2837 Mat B; 2838 PetscErrorCode ierr; 2839 PetscInt i,nz,header[4],*rowlengths = 0,M,N; 2840 int fd; 2841 PetscMPIInt size; 2842 MPI_Comm comm; 2843 2844 PetscFunctionBegin; 2845 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2846 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2847 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2848 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2849 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2850 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2851 M = header[1]; N = header[2]; nz = header[3]; 2852 2853 if (nz < 0) { 2854 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2855 } 2856 2857 /* read in row lengths */ 2858 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2859 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2860 2861 /* create our matrix */ 2862 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2863 ierr = MatSetType(B,type);CHKERRQ(ierr); 2864 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2865 a = (Mat_SeqAIJ*)B->data; 2866 2867 /* read in column indices and adjust for Fortran indexing*/ 2868 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2869 2870 /* read in nonzero values */ 2871 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2872 2873 /* set matrix "i" values */ 2874 a->i[0] = 0; 2875 for (i=1; i<= M; i++) { 2876 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2877 a->ilen[i-1] = rowlengths[i-1]; 2878 } 2879 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2880 2881 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2882 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2883 *A = B; 2884 PetscFunctionReturn(0); 2885 } 2886 2887 #undef __FUNCT__ 2888 #define __FUNCT__ "MatEqual_SeqAIJ" 2889 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2890 { 2891 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2892 PetscErrorCode ierr; 2893 2894 PetscFunctionBegin; 2895 /* If the matrix dimensions are not equal,or no of nonzeros */ 2896 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2897 *flg = PETSC_FALSE; 2898 PetscFunctionReturn(0); 2899 } 2900 2901 /* if the a->i are the same */ 2902 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2903 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2904 2905 /* if a->j are the same */ 2906 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2907 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2908 2909 /* if a->a are the same */ 2910 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2911 2912 PetscFunctionReturn(0); 2913 2914 } 2915 2916 #undef __FUNCT__ 2917 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2918 /*@C 2919 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2920 provided by the user. 2921 2922 Coolective on MPI_Comm 2923 2924 Input Parameters: 2925 + comm - must be an MPI communicator of size 1 2926 . m - number of rows 2927 . n - number of columns 2928 . i - row indices 2929 . j - column indices 2930 - a - matrix values 2931 2932 Output Parameter: 2933 . mat - the matrix 2934 2935 Level: intermediate 2936 2937 Notes: 2938 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2939 once the matrix is destroyed 2940 2941 You cannot set new nonzero locations into this matrix, that will generate an error. 2942 2943 The i and j indices are 0 based 2944 2945 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2946 2947 @*/ 2948 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2949 { 2950 PetscErrorCode ierr; 2951 PetscInt ii; 2952 Mat_SeqAIJ *aij; 2953 2954 PetscFunctionBegin; 2955 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 2956 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 2957 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);CHKERRQ(ierr); 2958 aij = (Mat_SeqAIJ*)(*mat)->data; 2959 2960 if (i[0] != 0) { 2961 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2962 } 2963 aij->i = i; 2964 aij->j = j; 2965 aij->a = a; 2966 aij->singlemalloc = PETSC_FALSE; 2967 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2968 aij->freedata = PETSC_FALSE; 2969 2970 for (ii=0; ii<m; ii++) { 2971 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2972 #if defined(PETSC_USE_BOPT_g) 2973 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2974 #endif 2975 } 2976 #if defined(PETSC_USE_BOPT_g) 2977 for (ii=0; ii<aij->i[m]; ii++) { 2978 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2979 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2980 } 2981 #endif 2982 2983 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2984 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2985 PetscFunctionReturn(0); 2986 } 2987 2988 #undef __FUNCT__ 2989 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2990 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2991 { 2992 PetscErrorCode ierr; 2993 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2994 2995 PetscFunctionBegin; 2996 if (coloring->ctype == IS_COLORING_LOCAL) { 2997 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2998 a->coloring = coloring; 2999 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3000 PetscInt i,*larray; 3001 ISColoring ocoloring; 3002 ISColoringValue *colors; 3003 3004 /* set coloring for diagonal portion */ 3005 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3006 for (i=0; i<A->n; i++) { 3007 larray[i] = i; 3008 } 3009 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3010 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3011 for (i=0; i<A->n; i++) { 3012 colors[i] = coloring->colors[larray[i]]; 3013 } 3014 ierr = PetscFree(larray);CHKERRQ(ierr); 3015 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3016 a->coloring = ocoloring; 3017 } 3018 PetscFunctionReturn(0); 3019 } 3020 3021 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3022 EXTERN_C_BEGIN 3023 #include "adic/ad_utils.h" 3024 EXTERN_C_END 3025 3026 #undef __FUNCT__ 3027 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3028 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3029 { 3030 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3031 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3032 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3033 ISColoringValue *color; 3034 3035 PetscFunctionBegin; 3036 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3037 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3038 color = a->coloring->colors; 3039 /* loop over rows */ 3040 for (i=0; i<m; i++) { 3041 nz = ii[i+1] - ii[i]; 3042 /* loop over columns putting computed value into matrix */ 3043 for (j=0; j<nz; j++) { 3044 *v++ = values[color[*jj++]]; 3045 } 3046 values += nlen; /* jump to next row of derivatives */ 3047 } 3048 PetscFunctionReturn(0); 3049 } 3050 3051 #else 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3055 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3056 { 3057 PetscFunctionBegin; 3058 SETERRQ(PETSC_ERR_SUP_SYS,"PETSc installed without ADIC"); 3059 } 3060 3061 #endif 3062 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3065 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3066 { 3067 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3068 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3069 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3070 ISColoringValue *color; 3071 3072 PetscFunctionBegin; 3073 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3074 color = a->coloring->colors; 3075 /* loop over rows */ 3076 for (i=0; i<m; i++) { 3077 nz = ii[i+1] - ii[i]; 3078 /* loop over columns putting computed value into matrix */ 3079 for (j=0; j<nz; j++) { 3080 *v++ = values[color[*jj++]]; 3081 } 3082 values += nl; /* jump to next row of derivatives */ 3083 } 3084 PetscFunctionReturn(0); 3085 } 3086 3087