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