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