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