1 #define PETSCMAT_DLL 2 3 4 /* 5 Defines the basic matrix operations for the AIJ (compressed row) 6 matrix storage format. 7 */ 8 9 10 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 11 #include "petscblaslapack.h" 12 #include "petscbt.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->rmap->n; 21 MatScalar *aa = aij->a; 22 PetscScalar *v; 23 PetscTruth missing; 24 25 PetscFunctionBegin; 26 if (Y->assembled) { 27 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 28 if (!missing) { 29 diag = aij->diag; 30 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 31 if (is == INSERT_VALUES) { 32 for (i=0; i<m; i++) { 33 aa[diag[i]] = v[i]; 34 } 35 } else { 36 for (i=0; i<m; i++) { 37 aa[diag[i]] += v[i]; 38 } 39 } 40 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 } 44 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 50 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 51 { 52 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 53 PetscErrorCode ierr; 54 PetscInt i,ishift; 55 56 PetscFunctionBegin; 57 *m = A->rmap->n; 58 if (!ia) PetscFunctionReturn(0); 59 ishift = 0; 60 if (symmetric && !A->structurally_symmetric) { 61 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 62 } else if (oshift == 1) { 63 PetscInt nz = a->i[A->rmap->n]; 64 /* malloc space and add 1 to i and j indices */ 65 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 66 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 67 if (ja) { 68 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 69 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 70 } 71 } else { 72 *ia = a->i; 73 if (ja) *ja = a->j; 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 80 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 if (!ia) PetscFunctionReturn(0); 86 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 87 ierr = PetscFree(*ia);CHKERRQ(ierr); 88 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 89 } 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 95 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 96 { 97 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 98 PetscErrorCode ierr; 99 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 100 PetscInt nz = a->i[m],row,*jj,mr,col; 101 102 PetscFunctionBegin; 103 *nn = n; 104 if (!ia) PetscFunctionReturn(0); 105 if (symmetric) { 106 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 107 } else { 108 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 109 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 110 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 111 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 112 jj = a->j; 113 for (i=0; i<nz; i++) { 114 collengths[jj[i]]++; 115 } 116 cia[0] = oshift; 117 for (i=0; i<n; i++) { 118 cia[i+1] = cia[i] + collengths[i]; 119 } 120 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 121 jj = a->j; 122 for (row=0; row<m; row++) { 123 mr = a->i[row+1] - a->i[row]; 124 for (i=0; i<mr; i++) { 125 col = *jj++; 126 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 127 } 128 } 129 ierr = PetscFree(collengths);CHKERRQ(ierr); 130 *ia = cia; *ja = cja; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 137 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 138 { 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (!ia) PetscFunctionReturn(0); 143 144 ierr = PetscFree(*ia);CHKERRQ(ierr); 145 ierr = PetscFree(*ja);CHKERRQ(ierr); 146 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 152 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 153 { 154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 155 PetscInt *ai = a->i; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatSetValues_SeqAIJ" 165 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 166 { 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 169 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 170 PetscErrorCode ierr; 171 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 172 MatScalar *ap,value,*aa = a->a; 173 PetscTruth ignorezeroentries = a->ignorezeroentries; 174 PetscTruth roworiented = a->roworiented; 175 176 PetscFunctionBegin; 177 if (v) PetscValidScalarPointer(v,6); 178 for (k=0; k<m; k++) { /* loop over added rows */ 179 row = im[k]; 180 if (row < 0) continue; 181 #if defined(PETSC_USE_DEBUG) 182 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 183 #endif 184 rp = aj + ai[row]; ap = aa + ai[row]; 185 rmax = imax[row]; nrow = ailen[row]; 186 low = 0; 187 high = nrow; 188 for (l=0; l<n; l++) { /* loop over added columns */ 189 if (in[l] < 0) continue; 190 #if defined(PETSC_USE_DEBUG) 191 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 192 #endif 193 col = in[l]; 194 if (v) { 195 if (roworiented) { 196 value = v[l + k*n]; 197 } else { 198 value = v[k + l*m]; 199 } 200 } else { 201 value = 0.; 202 } 203 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 204 205 if (col <= lastcol) low = 0; else high = nrow; 206 lastcol = col; 207 while (high-low > 5) { 208 t = (low+high)/2; 209 if (rp[t] > col) high = t; 210 else low = t; 211 } 212 for (i=low; i<high; i++) { 213 if (rp[i] > col) break; 214 if (rp[i] == col) { 215 if (is == ADD_VALUES) ap[i] += value; 216 else ap[i] = value; 217 low = i + 1; 218 goto noinsert; 219 } 220 } 221 if (value == 0.0 && ignorezeroentries) goto noinsert; 222 if (nonew == 1) goto noinsert; 223 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 224 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 225 N = nrow++ - 1; a->nz++; high++; 226 /* shift up all the later entries in this row */ 227 for (ii=N; ii>=i; ii--) { 228 rp[ii+1] = rp[ii]; 229 ap[ii+1] = ap[ii]; 230 } 231 rp[i] = col; 232 ap[i] = value; 233 low = i + 1; 234 noinsert:; 235 } 236 ailen[row] = nrow; 237 } 238 A->same_nonzero = PETSC_FALSE; 239 PetscFunctionReturn(0); 240 } 241 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatGetValues_SeqAIJ" 245 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 246 { 247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 248 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 249 PetscInt *ai = a->i,*ailen = a->ilen; 250 MatScalar *ap,*aa = a->a; 251 252 PetscFunctionBegin; 253 for (k=0; k<m; k++) { /* loop over rows */ 254 row = im[k]; 255 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 256 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 257 rp = aj + ai[row]; ap = aa + ai[row]; 258 nrow = ailen[row]; 259 for (l=0; l<n; l++) { /* loop over columns */ 260 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 261 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 262 col = in[l] ; 263 high = nrow; low = 0; /* assume unsorted */ 264 while (high-low > 5) { 265 t = (low+high)/2; 266 if (rp[t] > col) high = t; 267 else low = t; 268 } 269 for (i=low; i<high; i++) { 270 if (rp[i] > col) break; 271 if (rp[i] == col) { 272 *v++ = ap[i]; 273 goto finished; 274 } 275 } 276 *v++ = 0.0; 277 finished:; 278 } 279 } 280 PetscFunctionReturn(0); 281 } 282 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "MatView_SeqAIJ_Binary" 286 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 287 { 288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 289 PetscErrorCode ierr; 290 PetscInt i,*col_lens; 291 int fd; 292 293 PetscFunctionBegin; 294 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 295 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 296 col_lens[0] = MAT_FILE_CLASSID; 297 col_lens[1] = A->rmap->n; 298 col_lens[2] = A->cmap->n; 299 col_lens[3] = a->nz; 300 301 /* store lengths of each row and write (including header) to file */ 302 for (i=0; i<A->rmap->n; i++) { 303 col_lens[4+i] = a->i[i+1] - a->i[i]; 304 } 305 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 306 ierr = PetscFree(col_lens);CHKERRQ(ierr); 307 308 /* store column indices (zero start index) */ 309 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 310 311 /* store nonzero values */ 312 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 320 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 321 { 322 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 323 PetscErrorCode ierr; 324 PetscInt i,j,m = A->rmap->n,shift=0; 325 const char *name; 326 PetscViewerFormat format; 327 328 PetscFunctionBegin; 329 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 330 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 331 if (format == PETSC_VIEWER_ASCII_MATLAB) { 332 PetscInt nofinalvalue = 0; 333 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 334 nofinalvalue = 1; 335 } 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 341 342 for (i=0; i<m; i++) { 343 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 344 #if defined(PETSC_USE_COMPLEX) 345 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); 346 #else 347 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 348 #endif 349 } 350 } 351 if (nofinalvalue) { 352 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 353 } 354 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 355 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 356 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 357 PetscFunctionReturn(0); 358 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 360 for (i=0; i<m; i++) { 361 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 362 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 363 #if defined(PETSC_USE_COMPLEX) 364 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 365 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 366 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 367 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 368 } else if (PetscRealPart(a->a[j]) != 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 370 } 371 #else 372 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 373 #endif 374 } 375 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 376 } 377 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 378 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 379 PetscInt nzd=0,fshift=1,*sptr; 380 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 381 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 382 for (i=0; i<m; i++) { 383 sptr[i] = nzd+1; 384 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 385 if (a->j[j] >= i) { 386 #if defined(PETSC_USE_COMPLEX) 387 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 388 #else 389 if (a->a[j] != 0.0) nzd++; 390 #endif 391 } 392 } 393 } 394 sptr[m] = nzd+1; 395 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 396 for (i=0; i<m+1; i+=6) { 397 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);} 398 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);} 399 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);} 400 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 401 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 402 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 403 } 404 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 405 ierr = PetscFree(sptr);CHKERRQ(ierr); 406 for (i=0; i<m; i++) { 407 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 408 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 413 for (i=0; i<m; i++) { 414 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 415 if (a->j[j] >= i) { 416 #if defined(PETSC_USE_COMPLEX) 417 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 418 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 419 } 420 #else 421 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 422 #endif 423 } 424 } 425 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 426 } 427 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 428 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 429 PetscInt cnt = 0,jcnt; 430 PetscScalar value; 431 432 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 433 for (i=0; i<m; i++) { 434 jcnt = 0; 435 for (j=0; j<A->cmap->n; j++) { 436 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 437 value = a->a[cnt++]; 438 jcnt++; 439 } else { 440 value = 0.0; 441 } 442 #if defined(PETSC_USE_COMPLEX) 443 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 444 #else 445 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 446 #endif 447 } 448 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 449 } 450 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 451 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 452 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 453 #if defined(PETSC_USE_COMPLEX) 454 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 457 #endif 458 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 459 for (i=0; i<m; i++) { 460 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 461 #if defined(PETSC_USE_COMPLEX) 462 if (PetscImaginaryPart(a->a[j]) > 0.0) { 463 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 464 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 465 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 466 } else { 467 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 468 } 469 #else 470 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 471 #endif 472 } 473 } 474 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 475 } else { 476 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 477 if (A->factortype){ 478 for (i=0; i<m; i++) { 479 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 480 /* L part */ 481 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 482 #if defined(PETSC_USE_COMPLEX) 483 if (PetscImaginaryPart(a->a[j]) > 0.0) { 484 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 485 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 486 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 487 } else { 488 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 489 } 490 #else 491 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 492 #endif 493 } 494 /* diagonal */ 495 j = a->diag[i]; 496 #if defined(PETSC_USE_COMPLEX) 497 if (PetscImaginaryPart(a->a[j]) > 0.0) { 498 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 499 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 500 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 501 } else { 502 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 503 } 504 #else 505 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 506 #endif 507 508 /* U part */ 509 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 510 #if defined(PETSC_USE_COMPLEX) 511 if (PetscImaginaryPart(a->a[j]) > 0.0) { 512 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 513 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 514 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 515 } else { 516 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 517 } 518 #else 519 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 520 #endif 521 } 522 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 523 } 524 } else { 525 for (i=0; i<m; i++) { 526 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 527 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 528 #if defined(PETSC_USE_COMPLEX) 529 if (PetscImaginaryPart(a->a[j]) > 0.0) { 530 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 531 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 532 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 533 } else { 534 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 535 } 536 #else 537 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 538 #endif 539 } 540 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 541 } 542 } 543 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 544 } 545 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 551 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 552 { 553 Mat A = (Mat) Aa; 554 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 555 PetscErrorCode ierr; 556 PetscInt i,j,m = A->rmap->n,color; 557 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 558 PetscViewer viewer; 559 PetscViewerFormat format; 560 561 PetscFunctionBegin; 562 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 563 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 564 565 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 566 /* loop over matrix elements drawing boxes */ 567 568 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 569 /* Blue for negative, Cyan for zero and Red for positive */ 570 color = PETSC_DRAW_BLUE; 571 for (i=0; i<m; i++) { 572 y_l = m - i - 1.0; y_r = y_l + 1.0; 573 for (j=a->i[i]; j<a->i[i+1]; j++) { 574 x_l = a->j[j] ; x_r = x_l + 1.0; 575 #if defined(PETSC_USE_COMPLEX) 576 if (PetscRealPart(a->a[j]) >= 0.) continue; 577 #else 578 if (a->a[j] >= 0.) continue; 579 #endif 580 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 581 } 582 } 583 color = PETSC_DRAW_CYAN; 584 for (i=0; i<m; i++) { 585 y_l = m - i - 1.0; y_r = y_l + 1.0; 586 for (j=a->i[i]; j<a->i[i+1]; j++) { 587 x_l = a->j[j]; x_r = x_l + 1.0; 588 if (a->a[j] != 0.) continue; 589 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 590 } 591 } 592 color = PETSC_DRAW_RED; 593 for (i=0; i<m; i++) { 594 y_l = m - i - 1.0; y_r = y_l + 1.0; 595 for (j=a->i[i]; j<a->i[i+1]; j++) { 596 x_l = a->j[j]; x_r = x_l + 1.0; 597 #if defined(PETSC_USE_COMPLEX) 598 if (PetscRealPart(a->a[j]) <= 0.) continue; 599 #else 600 if (a->a[j] <= 0.) continue; 601 #endif 602 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 603 } 604 } 605 } else { 606 /* use contour shading to indicate magnitude of values */ 607 /* first determine max of all nonzero values */ 608 PetscInt nz = a->nz,count; 609 PetscDraw popup; 610 PetscReal scale; 611 612 for (i=0; i<nz; i++) { 613 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 614 } 615 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 616 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 617 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 618 count = 0; 619 for (i=0; i<m; i++) { 620 y_l = m - i - 1.0; y_r = y_l + 1.0; 621 for (j=a->i[i]; j<a->i[i+1]; j++) { 622 x_l = a->j[j]; x_r = x_l + 1.0; 623 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 624 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 625 count++; 626 } 627 } 628 } 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "MatView_SeqAIJ_Draw" 634 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 635 { 636 PetscErrorCode ierr; 637 PetscDraw draw; 638 PetscReal xr,yr,xl,yl,h,w; 639 PetscTruth isnull; 640 641 PetscFunctionBegin; 642 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 643 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 644 if (isnull) PetscFunctionReturn(0); 645 646 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 647 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 648 xr += w; yr += h; xl = -w; yl = -h; 649 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 650 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 651 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "MatView_SeqAIJ" 657 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 658 { 659 PetscErrorCode ierr; 660 PetscTruth iascii,isbinary,isdraw; 661 662 PetscFunctionBegin; 663 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 664 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 665 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 666 if (iascii) { 667 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 668 } else if (isbinary) { 669 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 670 } else if (isdraw) { 671 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 672 } else { 673 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 674 } 675 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 681 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 682 { 683 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 684 PetscErrorCode ierr; 685 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 686 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 687 MatScalar *aa = a->a,*ap; 688 PetscReal ratio=0.6; 689 690 PetscFunctionBegin; 691 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 692 693 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 694 for (i=1; i<m; i++) { 695 /* move each row back by the amount of empty slots (fshift) before it*/ 696 fshift += imax[i-1] - ailen[i-1]; 697 rmax = PetscMax(rmax,ailen[i]); 698 if (fshift) { 699 ip = aj + ai[i] ; 700 ap = aa + ai[i] ; 701 N = ailen[i]; 702 for (j=0; j<N; j++) { 703 ip[j-fshift] = ip[j]; 704 ap[j-fshift] = ap[j]; 705 } 706 } 707 ai[i] = ai[i-1] + ailen[i-1]; 708 } 709 if (m) { 710 fshift += imax[m-1] - ailen[m-1]; 711 ai[m] = ai[m-1] + ailen[m-1]; 712 } 713 /* reset ilen and imax for each row */ 714 for (i=0; i<m; i++) { 715 ailen[i] = imax[i] = ai[i+1] - ai[i]; 716 } 717 a->nz = ai[m]; 718 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 719 720 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 721 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 722 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 723 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 724 A->info.mallocs += a->reallocs; 725 a->reallocs = 0; 726 A->info.nz_unneeded = (double)fshift; 727 a->rmax = rmax; 728 729 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 730 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 731 A->same_nonzero = PETSC_TRUE; 732 733 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 734 735 a->idiagvalid = PETSC_FALSE; 736 PetscFunctionReturn(0); 737 } 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "MatRealPart_SeqAIJ" 741 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 742 { 743 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 744 PetscInt i,nz = a->nz; 745 MatScalar *aa = a->a; 746 747 PetscFunctionBegin; 748 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 754 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 755 { 756 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 757 PetscInt i,nz = a->nz; 758 MatScalar *aa = a->a; 759 760 PetscFunctionBegin; 761 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 762 PetscFunctionReturn(0); 763 } 764 765 #undef __FUNCT__ 766 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 767 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 768 { 769 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 774 PetscFunctionReturn(0); 775 } 776 777 #undef __FUNCT__ 778 #define __FUNCT__ "MatDestroy_SeqAIJ" 779 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 780 { 781 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 782 PetscErrorCode ierr; 783 784 PetscFunctionBegin; 785 #if defined(PETSC_USE_LOG) 786 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 787 #endif 788 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 789 if (a->row) { 790 ierr = ISDestroy(a->row);CHKERRQ(ierr); 791 } 792 if (a->col) { 793 ierr = ISDestroy(a->col);CHKERRQ(ierr); 794 } 795 ierr = PetscFree(a->diag);CHKERRQ(ierr); 796 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 797 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 798 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 799 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 800 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 801 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 802 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 803 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 804 if (a->compressedrow.checked && a->compressedrow.use){ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);} 805 806 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 807 808 ierr = PetscFree(a);CHKERRQ(ierr); 809 810 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 811 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 812 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 813 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 814 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 815 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 817 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "MatSetOption_SeqAIJ" 826 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 827 { 828 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 switch (op) { 833 case MAT_ROW_ORIENTED: 834 a->roworiented = flg; 835 break; 836 case MAT_KEEP_NONZERO_PATTERN: 837 a->keepnonzeropattern = flg; 838 break; 839 case MAT_NEW_NONZERO_LOCATIONS: 840 a->nonew = (flg ? 0 : 1); 841 break; 842 case MAT_NEW_NONZERO_LOCATION_ERR: 843 a->nonew = (flg ? -1 : 0); 844 break; 845 case MAT_NEW_NONZERO_ALLOCATION_ERR: 846 a->nonew = (flg ? -2 : 0); 847 break; 848 case MAT_UNUSED_NONZERO_LOCATION_ERR: 849 a->nounused = (flg ? -1 : 0); 850 break; 851 case MAT_IGNORE_ZERO_ENTRIES: 852 a->ignorezeroentries = flg; 853 break; 854 case MAT_USE_COMPRESSEDROW: 855 a->compressedrow.use = flg; 856 break; 857 case MAT_SPD: 858 A->spd_set = PETSC_TRUE; 859 A->spd = flg; 860 if (flg) { 861 A->symmetric = PETSC_TRUE; 862 A->structurally_symmetric = PETSC_TRUE; 863 A->symmetric_set = PETSC_TRUE; 864 A->structurally_symmetric_set = PETSC_TRUE; 865 } 866 break; 867 case MAT_SYMMETRIC: 868 case MAT_STRUCTURALLY_SYMMETRIC: 869 case MAT_HERMITIAN: 870 case MAT_SYMMETRY_ETERNAL: 871 case MAT_NEW_DIAGONALS: 872 case MAT_IGNORE_OFF_PROC_ENTRIES: 873 case MAT_USE_HASH_TABLE: 874 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 875 break; 876 case MAT_USE_INODES: 877 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 878 break; 879 default: 880 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 881 } 882 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 #undef __FUNCT__ 887 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 888 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 889 { 890 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 891 PetscErrorCode ierr; 892 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 893 PetscScalar *aa=a->a,*x,zero=0.0; 894 895 PetscFunctionBegin; 896 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 897 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 898 899 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU){ 900 PetscInt *diag=a->diag; 901 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 902 for (i=0; i<n; i++) x[i] = aa[diag[i]]; 903 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 904 PetscFunctionReturn(0); 905 } 906 907 ierr = VecSet(v,zero);CHKERRQ(ierr); 908 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 909 for (i=0; i<n; i++) { 910 nz = ai[i+1] - ai[i]; 911 if (!nz) x[i] = 0.0; 912 for (j=ai[i]; j<ai[i+1]; j++){ 913 if (aj[j] == i) { 914 x[i] = aa[j]; 915 break; 916 } 917 } 918 } 919 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 920 PetscFunctionReturn(0); 921 } 922 923 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 924 #undef __FUNCT__ 925 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 926 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 927 { 928 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 929 PetscScalar *x,*y; 930 PetscErrorCode ierr; 931 PetscInt m = A->rmap->n; 932 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 933 MatScalar *v; 934 PetscScalar alpha; 935 PetscInt n,i,j,*idx,*ii,*ridx=PETSC_NULL; 936 Mat_CompressedRow cprow = a->compressedrow; 937 PetscTruth usecprow = cprow.use; 938 #endif 939 940 PetscFunctionBegin; 941 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 942 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 943 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 944 945 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 946 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 947 #else 948 if (usecprow){ 949 m = cprow.nrows; 950 ii = cprow.i; 951 ridx = cprow.rindex; 952 } else { 953 ii = a->i; 954 } 955 for (i=0; i<m; i++) { 956 idx = a->j + ii[i] ; 957 v = a->a + ii[i] ; 958 n = ii[i+1] - ii[i]; 959 if (usecprow){ 960 alpha = x[ridx[i]]; 961 } else { 962 alpha = x[i]; 963 } 964 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 965 } 966 #endif 967 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 968 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 969 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 970 PetscFunctionReturn(0); 971 } 972 973 #undef __FUNCT__ 974 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 975 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 976 { 977 PetscErrorCode ierr; 978 979 PetscFunctionBegin; 980 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 981 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 982 PetscFunctionReturn(0); 983 } 984 985 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 986 #undef __FUNCT__ 987 #define __FUNCT__ "MatMult_SeqAIJ" 988 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 989 { 990 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 991 PetscScalar *y; 992 const PetscScalar *x; 993 const MatScalar *aa; 994 PetscErrorCode ierr; 995 PetscInt m=A->rmap->n; 996 const PetscInt *aj,*ii,*ridx=PETSC_NULL; 997 PetscInt n,i,nonzerorow=0; 998 PetscScalar sum; 999 PetscTruth usecprow=a->compressedrow.use; 1000 1001 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1002 #pragma disjoint(*x,*y,*aa) 1003 #endif 1004 1005 PetscFunctionBegin; 1006 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1007 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1008 aj = a->j; 1009 aa = a->a; 1010 ii = a->i; 1011 if (usecprow){ /* use compressed row format */ 1012 m = a->compressedrow.nrows; 1013 ii = a->compressedrow.i; 1014 ridx = a->compressedrow.rindex; 1015 for (i=0; i<m; i++){ 1016 n = ii[i+1] - ii[i]; 1017 aj = a->j + ii[i]; 1018 aa = a->a + ii[i]; 1019 sum = 0.0; 1020 nonzerorow += (n>0); 1021 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1022 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1023 y[*ridx++] = sum; 1024 } 1025 } else { /* do not use compressed row format */ 1026 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1027 fortranmultaij_(&m,x,ii,aj,aa,y); 1028 #else 1029 for (i=0; i<m; i++) { 1030 n = ii[i+1] - ii[i]; 1031 aj = a->j + ii[i]; 1032 aa = a->a + ii[i]; 1033 sum = 0.0; 1034 nonzerorow += (n>0); 1035 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1036 y[i] = sum; 1037 } 1038 #endif 1039 } 1040 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1041 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1042 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1043 PetscFunctionReturn(0); 1044 } 1045 1046 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1049 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1050 { 1051 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1052 PetscScalar *x,*y,*z; 1053 const MatScalar *aa; 1054 PetscErrorCode ierr; 1055 PetscInt m = A->rmap->n,*aj,*ii; 1056 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1057 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1058 PetscScalar sum; 1059 PetscTruth usecprow=a->compressedrow.use; 1060 #endif 1061 1062 PetscFunctionBegin; 1063 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1064 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1065 if (zz != yy) { 1066 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1067 } else { 1068 z = y; 1069 } 1070 1071 aj = a->j; 1072 aa = a->a; 1073 ii = a->i; 1074 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1075 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1076 #else 1077 if (usecprow){ /* use compressed row format */ 1078 if (zz != yy){ 1079 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1080 } 1081 m = a->compressedrow.nrows; 1082 ii = a->compressedrow.i; 1083 ridx = a->compressedrow.rindex; 1084 for (i=0; i<m; i++){ 1085 n = ii[i+1] - ii[i]; 1086 aj = a->j + ii[i]; 1087 aa = a->a + ii[i]; 1088 sum = y[*ridx]; 1089 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1090 z[*ridx++] = sum; 1091 } 1092 } else { /* do not use compressed row format */ 1093 for (i=0; i<m; i++) { 1094 jrow = ii[i]; 1095 n = ii[i+1] - jrow; 1096 sum = y[i]; 1097 for (j=0; j<n; j++) { 1098 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1099 } 1100 z[i] = sum; 1101 } 1102 } 1103 #endif 1104 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1105 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1106 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1107 if (zz != yy) { 1108 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1109 } 1110 PetscFunctionReturn(0); 1111 } 1112 1113 /* 1114 Adds diagonal pointers to sparse matrix structure. 1115 */ 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1118 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1119 { 1120 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1121 PetscErrorCode ierr; 1122 PetscInt i,j,m = A->rmap->n; 1123 1124 PetscFunctionBegin; 1125 if (!a->diag) { 1126 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1127 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1128 } 1129 for (i=0; i<A->rmap->n; i++) { 1130 a->diag[i] = a->i[i+1]; 1131 for (j=a->i[i]; j<a->i[i+1]; j++) { 1132 if (a->j[j] == i) { 1133 a->diag[i] = j; 1134 break; 1135 } 1136 } 1137 } 1138 PetscFunctionReturn(0); 1139 } 1140 1141 /* 1142 Checks for missing diagonals 1143 */ 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1146 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1147 { 1148 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1149 PetscInt *diag,*jj = a->j,i; 1150 1151 PetscFunctionBegin; 1152 *missing = PETSC_FALSE; 1153 if (A->rmap->n > 0 && !jj) { 1154 *missing = PETSC_TRUE; 1155 if (d) *d = 0; 1156 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1157 } else { 1158 diag = a->diag; 1159 for (i=0; i<A->rmap->n; i++) { 1160 if (jj[diag[i]] != i) { 1161 *missing = PETSC_TRUE; 1162 if (d) *d = i; 1163 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1164 } 1165 } 1166 } 1167 PetscFunctionReturn(0); 1168 } 1169 1170 EXTERN_C_BEGIN 1171 #undef __FUNCT__ 1172 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1173 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1174 { 1175 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1176 PetscErrorCode ierr; 1177 PetscInt i,*diag,m = A->rmap->n; 1178 MatScalar *v = a->a; 1179 PetscScalar *idiag,*mdiag; 1180 1181 PetscFunctionBegin; 1182 if (a->idiagvalid) PetscFunctionReturn(0); 1183 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1184 diag = a->diag; 1185 if (!a->idiag) { 1186 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1187 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1188 v = a->a; 1189 } 1190 mdiag = a->mdiag; 1191 idiag = a->idiag; 1192 1193 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1194 for (i=0; i<m; i++) { 1195 mdiag[i] = v[diag[i]]; 1196 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1197 idiag[i] = 1.0/v[diag[i]]; 1198 } 1199 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1200 } else { 1201 for (i=0; i<m; i++) { 1202 mdiag[i] = v[diag[i]]; 1203 idiag[i] = omega/(fshift + v[diag[i]]); 1204 } 1205 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1206 } 1207 a->idiagvalid = PETSC_TRUE; 1208 PetscFunctionReturn(0); 1209 } 1210 EXTERN_C_END 1211 1212 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "MatSOR_SeqAIJ" 1215 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1216 { 1217 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1218 PetscScalar *x,d,sum,*t,scale; 1219 const MatScalar *v = a->a,*idiag=0,*mdiag; 1220 const PetscScalar *b, *bs,*xb, *ts; 1221 PetscErrorCode ierr; 1222 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1223 const PetscInt *idx,*diag; 1224 1225 PetscFunctionBegin; 1226 its = its*lits; 1227 1228 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1229 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1230 a->fshift = fshift; 1231 a->omega = omega; 1232 1233 diag = a->diag; 1234 t = a->ssor_work; 1235 idiag = a->idiag; 1236 mdiag = a->mdiag; 1237 1238 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1239 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1240 CHKMEMQ; 1241 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1242 if (flag == SOR_APPLY_UPPER) { 1243 /* apply (U + D/omega) to the vector */ 1244 bs = b; 1245 for (i=0; i<m; i++) { 1246 d = fshift + mdiag[i]; 1247 n = a->i[i+1] - diag[i] - 1; 1248 idx = a->j + diag[i] + 1; 1249 v = a->a + diag[i] + 1; 1250 sum = b[i]*d/omega; 1251 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1252 x[i] = sum; 1253 } 1254 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1255 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1256 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1257 PetscFunctionReturn(0); 1258 } 1259 1260 if (flag == SOR_APPLY_LOWER) { 1261 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1262 } else if (flag & SOR_EISENSTAT) { 1263 /* Let A = L + U + D; where L is lower trianglar, 1264 U is upper triangular, E = D/omega; This routine applies 1265 1266 (L + E)^{-1} A (U + E)^{-1} 1267 1268 to a vector efficiently using Eisenstat's trick. 1269 */ 1270 scale = (2.0/omega) - 1.0; 1271 1272 /* x = (E + U)^{-1} b */ 1273 for (i=m-1; i>=0; i--) { 1274 n = a->i[i+1] - diag[i] - 1; 1275 idx = a->j + diag[i] + 1; 1276 v = a->a + diag[i] + 1; 1277 sum = b[i]; 1278 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1279 x[i] = sum*idiag[i]; 1280 } 1281 1282 /* t = b - (2*E - D)x */ 1283 v = a->a; 1284 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1285 1286 /* t = (E + L)^{-1}t */ 1287 ts = t; 1288 diag = a->diag; 1289 for (i=0; i<m; i++) { 1290 n = diag[i] - a->i[i]; 1291 idx = a->j + a->i[i]; 1292 v = a->a + a->i[i]; 1293 sum = t[i]; 1294 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1295 t[i] = sum*idiag[i]; 1296 /* x = x + t */ 1297 x[i] += t[i]; 1298 } 1299 1300 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1301 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1302 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1303 PetscFunctionReturn(0); 1304 } 1305 if (flag & SOR_ZERO_INITIAL_GUESS) { 1306 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1307 for (i=0; i<m; i++) { 1308 n = diag[i] - a->i[i]; 1309 idx = a->j + a->i[i]; 1310 v = a->a + a->i[i]; 1311 sum = b[i]; 1312 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1313 t[i] = sum; 1314 x[i] = sum*idiag[i]; 1315 } 1316 xb = t; 1317 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1318 } else xb = b; 1319 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1320 for (i=m-1; i>=0; i--) { 1321 n = a->i[i+1] - diag[i] - 1; 1322 idx = a->j + diag[i] + 1; 1323 v = a->a + diag[i] + 1; 1324 sum = xb[i]; 1325 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1326 if (xb == b) { 1327 x[i] = sum*idiag[i]; 1328 } else { 1329 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1330 } 1331 } 1332 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1333 } 1334 its--; 1335 } 1336 while (its--) { 1337 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1338 for (i=0; i<m; i++) { 1339 n = a->i[i+1] - a->i[i]; 1340 idx = a->j + a->i[i]; 1341 v = a->a + a->i[i]; 1342 sum = b[i]; 1343 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1344 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1345 } 1346 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1347 } 1348 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1349 for (i=m-1; i>=0; i--) { 1350 n = a->i[i+1] - a->i[i]; 1351 idx = a->j + a->i[i]; 1352 v = a->a + a->i[i]; 1353 sum = b[i]; 1354 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1355 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1356 } 1357 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1358 } 1359 } 1360 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1361 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1362 CHKMEMQ; PetscFunctionReturn(0); 1363 } 1364 1365 1366 #undef __FUNCT__ 1367 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1368 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1369 { 1370 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1371 1372 PetscFunctionBegin; 1373 info->block_size = 1.0; 1374 info->nz_allocated = (double)a->maxnz; 1375 info->nz_used = (double)a->nz; 1376 info->nz_unneeded = (double)(a->maxnz - a->nz); 1377 info->assemblies = (double)A->num_ass; 1378 info->mallocs = (double)A->info.mallocs; 1379 info->memory = ((PetscObject)A)->mem; 1380 if (A->factortype) { 1381 info->fill_ratio_given = A->info.fill_ratio_given; 1382 info->fill_ratio_needed = A->info.fill_ratio_needed; 1383 info->factor_mallocs = A->info.factor_mallocs; 1384 } else { 1385 info->fill_ratio_given = 0; 1386 info->fill_ratio_needed = 0; 1387 info->factor_mallocs = 0; 1388 } 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1394 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1395 { 1396 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1397 PetscInt i,m = A->rmap->n - 1,d = 0; 1398 PetscErrorCode ierr; 1399 PetscTruth missing; 1400 1401 PetscFunctionBegin; 1402 if (a->keepnonzeropattern) { 1403 for (i=0; i<N; i++) { 1404 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1405 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1406 } 1407 if (diag != 0.0) { 1408 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1409 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1410 for (i=0; i<N; i++) { 1411 a->a[a->diag[rows[i]]] = diag; 1412 } 1413 } 1414 A->same_nonzero = PETSC_TRUE; 1415 } else { 1416 if (diag != 0.0) { 1417 for (i=0; i<N; i++) { 1418 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1419 if (a->ilen[rows[i]] > 0) { 1420 a->ilen[rows[i]] = 1; 1421 a->a[a->i[rows[i]]] = diag; 1422 a->j[a->i[rows[i]]] = rows[i]; 1423 } else { /* in case row was completely empty */ 1424 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1425 } 1426 } 1427 } else { 1428 for (i=0; i<N; i++) { 1429 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1430 a->ilen[rows[i]] = 0; 1431 } 1432 } 1433 A->same_nonzero = PETSC_FALSE; 1434 } 1435 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "MatGetRow_SeqAIJ" 1441 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1442 { 1443 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1444 PetscInt *itmp; 1445 1446 PetscFunctionBegin; 1447 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1448 1449 *nz = a->i[row+1] - a->i[row]; 1450 if (v) *v = a->a + a->i[row]; 1451 if (idx) { 1452 itmp = a->j + a->i[row]; 1453 if (*nz) { 1454 *idx = itmp; 1455 } 1456 else *idx = 0; 1457 } 1458 PetscFunctionReturn(0); 1459 } 1460 1461 /* remove this function? */ 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1464 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1465 { 1466 PetscFunctionBegin; 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "MatNorm_SeqAIJ" 1472 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1473 { 1474 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1475 MatScalar *v = a->a; 1476 PetscReal sum = 0.0; 1477 PetscErrorCode ierr; 1478 PetscInt i,j; 1479 1480 PetscFunctionBegin; 1481 if (type == NORM_FROBENIUS) { 1482 for (i=0; i<a->nz; i++) { 1483 #if defined(PETSC_USE_COMPLEX) 1484 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1485 #else 1486 sum += (*v)*(*v); v++; 1487 #endif 1488 } 1489 *nrm = sqrt(sum); 1490 } else if (type == NORM_1) { 1491 PetscReal *tmp; 1492 PetscInt *jj = a->j; 1493 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1494 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1495 *nrm = 0.0; 1496 for (j=0; j<a->nz; j++) { 1497 tmp[*jj++] += PetscAbsScalar(*v); v++; 1498 } 1499 for (j=0; j<A->cmap->n; j++) { 1500 if (tmp[j] > *nrm) *nrm = tmp[j]; 1501 } 1502 ierr = PetscFree(tmp);CHKERRQ(ierr); 1503 } else if (type == NORM_INFINITY) { 1504 *nrm = 0.0; 1505 for (j=0; j<A->rmap->n; j++) { 1506 v = a->a + a->i[j]; 1507 sum = 0.0; 1508 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1509 sum += PetscAbsScalar(*v); v++; 1510 } 1511 if (sum > *nrm) *nrm = sum; 1512 } 1513 } else { 1514 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1515 } 1516 PetscFunctionReturn(0); 1517 } 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "MatTranspose_SeqAIJ" 1521 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1522 { 1523 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1524 Mat C; 1525 PetscErrorCode ierr; 1526 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1527 MatScalar *array = a->a; 1528 1529 PetscFunctionBegin; 1530 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1531 1532 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1533 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1534 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1535 1536 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1537 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1538 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1539 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1540 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1541 ierr = PetscFree(col);CHKERRQ(ierr); 1542 } else { 1543 C = *B; 1544 } 1545 1546 for (i=0; i<m; i++) { 1547 len = ai[i+1]-ai[i]; 1548 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1549 array += len; 1550 aj += len; 1551 } 1552 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1553 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1554 1555 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1556 *B = C; 1557 } else { 1558 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1559 } 1560 PetscFunctionReturn(0); 1561 } 1562 1563 EXTERN_C_BEGIN 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1566 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1567 { 1568 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1569 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1570 MatScalar *va,*vb; 1571 PetscErrorCode ierr; 1572 PetscInt ma,na,mb,nb, i; 1573 1574 PetscFunctionBegin; 1575 bij = (Mat_SeqAIJ *) B->data; 1576 1577 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1578 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1579 if (ma!=nb || na!=mb){ 1580 *f = PETSC_FALSE; 1581 PetscFunctionReturn(0); 1582 } 1583 aii = aij->i; bii = bij->i; 1584 adx = aij->j; bdx = bij->j; 1585 va = aij->a; vb = bij->a; 1586 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1587 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1588 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1589 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1590 1591 *f = PETSC_TRUE; 1592 for (i=0; i<ma; i++) { 1593 while (aptr[i]<aii[i+1]) { 1594 PetscInt idc,idr; 1595 PetscScalar vc,vr; 1596 /* column/row index/value */ 1597 idc = adx[aptr[i]]; 1598 idr = bdx[bptr[idc]]; 1599 vc = va[aptr[i]]; 1600 vr = vb[bptr[idc]]; 1601 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1602 *f = PETSC_FALSE; 1603 goto done; 1604 } else { 1605 aptr[i]++; 1606 if (B || i!=idc) bptr[idc]++; 1607 } 1608 } 1609 } 1610 done: 1611 ierr = PetscFree(aptr);CHKERRQ(ierr); 1612 if (B) { 1613 ierr = PetscFree(bptr);CHKERRQ(ierr); 1614 } 1615 PetscFunctionReturn(0); 1616 } 1617 EXTERN_C_END 1618 1619 EXTERN_C_BEGIN 1620 #undef __FUNCT__ 1621 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1622 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1623 { 1624 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1625 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1626 MatScalar *va,*vb; 1627 PetscErrorCode ierr; 1628 PetscInt ma,na,mb,nb, i; 1629 1630 PetscFunctionBegin; 1631 bij = (Mat_SeqAIJ *) B->data; 1632 1633 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1634 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1635 if (ma!=nb || na!=mb){ 1636 *f = PETSC_FALSE; 1637 PetscFunctionReturn(0); 1638 } 1639 aii = aij->i; bii = bij->i; 1640 adx = aij->j; bdx = bij->j; 1641 va = aij->a; vb = bij->a; 1642 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1643 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1644 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1645 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1646 1647 *f = PETSC_TRUE; 1648 for (i=0; i<ma; i++) { 1649 while (aptr[i]<aii[i+1]) { 1650 PetscInt idc,idr; 1651 PetscScalar vc,vr; 1652 /* column/row index/value */ 1653 idc = adx[aptr[i]]; 1654 idr = bdx[bptr[idc]]; 1655 vc = va[aptr[i]]; 1656 vr = vb[bptr[idc]]; 1657 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1658 *f = PETSC_FALSE; 1659 goto done; 1660 } else { 1661 aptr[i]++; 1662 if (B || i!=idc) bptr[idc]++; 1663 } 1664 } 1665 } 1666 done: 1667 ierr = PetscFree(aptr);CHKERRQ(ierr); 1668 if (B) { 1669 ierr = PetscFree(bptr);CHKERRQ(ierr); 1670 } 1671 PetscFunctionReturn(0); 1672 } 1673 EXTERN_C_END 1674 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1677 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1678 { 1679 PetscErrorCode ierr; 1680 PetscFunctionBegin; 1681 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1682 PetscFunctionReturn(0); 1683 } 1684 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1687 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1688 { 1689 PetscErrorCode ierr; 1690 PetscFunctionBegin; 1691 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1697 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1698 { 1699 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1700 PetscScalar *l,*r,x; 1701 MatScalar *v; 1702 PetscErrorCode ierr; 1703 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1704 1705 PetscFunctionBegin; 1706 if (ll) { 1707 /* The local size is used so that VecMPI can be passed to this routine 1708 by MatDiagonalScale_MPIAIJ */ 1709 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1710 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1711 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1712 v = a->a; 1713 for (i=0; i<m; i++) { 1714 x = l[i]; 1715 M = a->i[i+1] - a->i[i]; 1716 for (j=0; j<M; j++) { (*v++) *= x;} 1717 } 1718 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1719 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1720 } 1721 if (rr) { 1722 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1723 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1724 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1725 v = a->a; jj = a->j; 1726 for (i=0; i<nz; i++) { 1727 (*v++) *= r[*jj++]; 1728 } 1729 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1730 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1731 } 1732 PetscFunctionReturn(0); 1733 } 1734 1735 #undef __FUNCT__ 1736 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1737 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1738 { 1739 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1740 PetscErrorCode ierr; 1741 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1742 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1743 const PetscInt *irow,*icol; 1744 PetscInt nrows,ncols; 1745 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1746 MatScalar *a_new,*mat_a; 1747 Mat C; 1748 PetscTruth stride,sorted; 1749 1750 PetscFunctionBegin; 1751 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1752 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1753 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1754 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1755 1756 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1757 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1758 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1759 1760 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1761 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1762 if (stride && step == 1) { 1763 /* special case of contiguous rows */ 1764 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1765 /* loop over new rows determining lens and starting points */ 1766 for (i=0; i<nrows; i++) { 1767 kstart = ai[irow[i]]; 1768 kend = kstart + ailen[irow[i]]; 1769 for (k=kstart; k<kend; k++) { 1770 if (aj[k] >= first) { 1771 starts[i] = k; 1772 break; 1773 } 1774 } 1775 sum = 0; 1776 while (k < kend) { 1777 if (aj[k++] >= first+ncols) break; 1778 sum++; 1779 } 1780 lens[i] = sum; 1781 } 1782 /* create submatrix */ 1783 if (scall == MAT_REUSE_MATRIX) { 1784 PetscInt n_cols,n_rows; 1785 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1786 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1787 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1788 C = *B; 1789 } else { 1790 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1791 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1792 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1793 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1794 } 1795 c = (Mat_SeqAIJ*)C->data; 1796 1797 /* loop over rows inserting into submatrix */ 1798 a_new = c->a; 1799 j_new = c->j; 1800 i_new = c->i; 1801 1802 for (i=0; i<nrows; i++) { 1803 ii = starts[i]; 1804 lensi = lens[i]; 1805 for (k=0; k<lensi; k++) { 1806 *j_new++ = aj[ii+k] - first; 1807 } 1808 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1809 a_new += lensi; 1810 i_new[i+1] = i_new[i] + lensi; 1811 c->ilen[i] = lensi; 1812 } 1813 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1814 } else { 1815 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1816 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1817 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1818 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1819 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1820 /* determine lens of each row */ 1821 for (i=0; i<nrows; i++) { 1822 kstart = ai[irow[i]]; 1823 kend = kstart + a->ilen[irow[i]]; 1824 lens[i] = 0; 1825 for (k=kstart; k<kend; k++) { 1826 if (smap[aj[k]]) { 1827 lens[i]++; 1828 } 1829 } 1830 } 1831 /* Create and fill new matrix */ 1832 if (scall == MAT_REUSE_MATRIX) { 1833 PetscTruth equal; 1834 1835 c = (Mat_SeqAIJ *)((*B)->data); 1836 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1837 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1838 if (!equal) { 1839 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1840 } 1841 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1842 C = *B; 1843 } else { 1844 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1845 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1846 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1847 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1848 } 1849 c = (Mat_SeqAIJ *)(C->data); 1850 for (i=0; i<nrows; i++) { 1851 row = irow[i]; 1852 kstart = ai[row]; 1853 kend = kstart + a->ilen[row]; 1854 mat_i = c->i[i]; 1855 mat_j = c->j + mat_i; 1856 mat_a = c->a + mat_i; 1857 mat_ilen = c->ilen + i; 1858 for (k=kstart; k<kend; k++) { 1859 if ((tcol=smap[a->j[k]])) { 1860 *mat_j++ = tcol - 1; 1861 *mat_a++ = a->a[k]; 1862 (*mat_ilen)++; 1863 1864 } 1865 } 1866 } 1867 /* Free work space */ 1868 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1869 ierr = PetscFree(smap);CHKERRQ(ierr); 1870 ierr = PetscFree(lens);CHKERRQ(ierr); 1871 } 1872 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1873 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1874 1875 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1876 *B = C; 1877 PetscFunctionReturn(0); 1878 } 1879 1880 #undef __FUNCT__ 1881 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1882 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1883 { 1884 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1885 PetscErrorCode ierr; 1886 Mat outA; 1887 PetscTruth row_identity,col_identity; 1888 1889 PetscFunctionBegin; 1890 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1891 1892 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1893 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1894 1895 outA = inA; 1896 outA->factortype = MAT_FACTOR_LU; 1897 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1898 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1899 a->row = row; 1900 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1901 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1902 a->col = col; 1903 1904 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1905 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1906 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1907 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1908 1909 if (!a->solve_work) { /* this matrix may have been factored before */ 1910 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1911 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1912 } 1913 1914 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1915 if (row_identity && col_identity) { 1916 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1917 } else { 1918 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1919 } 1920 PetscFunctionReturn(0); 1921 } 1922 1923 #undef __FUNCT__ 1924 #define __FUNCT__ "MatScale_SeqAIJ" 1925 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1926 { 1927 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1928 PetscScalar oalpha = alpha; 1929 PetscErrorCode ierr; 1930 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1931 1932 PetscFunctionBegin; 1933 BLASscal_(&bnz,&oalpha,a->a,&one); 1934 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1935 PetscFunctionReturn(0); 1936 } 1937 1938 #undef __FUNCT__ 1939 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1940 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1941 { 1942 PetscErrorCode ierr; 1943 PetscInt i; 1944 1945 PetscFunctionBegin; 1946 if (scall == MAT_INITIAL_MATRIX) { 1947 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1948 } 1949 1950 for (i=0; i<n; i++) { 1951 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1952 } 1953 PetscFunctionReturn(0); 1954 } 1955 1956 #undef __FUNCT__ 1957 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1958 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1959 { 1960 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1961 PetscErrorCode ierr; 1962 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1963 const PetscInt *idx; 1964 PetscInt start,end,*ai,*aj; 1965 PetscBT table; 1966 1967 PetscFunctionBegin; 1968 m = A->rmap->n; 1969 ai = a->i; 1970 aj = a->j; 1971 1972 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1973 1974 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1975 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1976 1977 for (i=0; i<is_max; i++) { 1978 /* Initialize the two local arrays */ 1979 isz = 0; 1980 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1981 1982 /* Extract the indices, assume there can be duplicate entries */ 1983 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1984 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1985 1986 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1987 for (j=0; j<n ; ++j){ 1988 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1989 } 1990 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1991 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1992 1993 k = 0; 1994 for (j=0; j<ov; j++){ /* for each overlap */ 1995 n = isz; 1996 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1997 row = nidx[k]; 1998 start = ai[row]; 1999 end = ai[row+1]; 2000 for (l = start; l<end ; l++){ 2001 val = aj[l] ; 2002 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2003 } 2004 } 2005 } 2006 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2007 } 2008 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2009 ierr = PetscFree(nidx);CHKERRQ(ierr); 2010 PetscFunctionReturn(0); 2011 } 2012 2013 /* -------------------------------------------------------------- */ 2014 #undef __FUNCT__ 2015 #define __FUNCT__ "MatPermute_SeqAIJ" 2016 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2017 { 2018 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2019 PetscErrorCode ierr; 2020 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2021 const PetscInt *row,*col; 2022 PetscInt *cnew,j,*lens; 2023 IS icolp,irowp; 2024 PetscInt *cwork = PETSC_NULL; 2025 PetscScalar *vwork = PETSC_NULL; 2026 2027 PetscFunctionBegin; 2028 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2029 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2030 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2031 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2032 2033 /* determine lengths of permuted rows */ 2034 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2035 for (i=0; i<m; i++) { 2036 lens[row[i]] = a->i[i+1] - a->i[i]; 2037 } 2038 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2039 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2040 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2041 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2042 ierr = PetscFree(lens);CHKERRQ(ierr); 2043 2044 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2045 for (i=0; i<m; i++) { 2046 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2047 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2048 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2049 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2050 } 2051 ierr = PetscFree(cnew);CHKERRQ(ierr); 2052 (*B)->assembled = PETSC_FALSE; 2053 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2054 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2055 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2056 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2057 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2058 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2059 PetscFunctionReturn(0); 2060 } 2061 2062 #undef __FUNCT__ 2063 #define __FUNCT__ "MatCopy_SeqAIJ" 2064 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2065 { 2066 PetscErrorCode ierr; 2067 2068 PetscFunctionBegin; 2069 /* If the two matrices have the same copy implementation, use fast copy. */ 2070 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2071 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2072 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2073 2074 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2075 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2076 } 2077 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2078 } else { 2079 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2080 } 2081 PetscFunctionReturn(0); 2082 } 2083 2084 #undef __FUNCT__ 2085 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2086 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2087 { 2088 PetscErrorCode ierr; 2089 2090 PetscFunctionBegin; 2091 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2092 PetscFunctionReturn(0); 2093 } 2094 2095 #undef __FUNCT__ 2096 #define __FUNCT__ "MatGetArray_SeqAIJ" 2097 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2098 { 2099 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2100 PetscFunctionBegin; 2101 *array = a->a; 2102 PetscFunctionReturn(0); 2103 } 2104 2105 #undef __FUNCT__ 2106 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2107 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2108 { 2109 PetscFunctionBegin; 2110 PetscFunctionReturn(0); 2111 } 2112 2113 #undef __FUNCT__ 2114 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2115 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2116 { 2117 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2118 PetscErrorCode ierr; 2119 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2120 PetscScalar dx,*y,*xx,*w3_array; 2121 PetscScalar *vscale_array; 2122 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2123 Vec w1,w2,w3; 2124 void *fctx = coloring->fctx; 2125 PetscTruth flg = PETSC_FALSE; 2126 2127 PetscFunctionBegin; 2128 if (!coloring->w1) { 2129 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2130 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2131 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2132 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2133 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2134 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2135 } 2136 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2137 2138 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2139 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2140 if (flg) { 2141 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2142 } else { 2143 PetscTruth assembled; 2144 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2145 if (assembled) { 2146 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2147 } 2148 } 2149 2150 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2151 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2152 2153 /* 2154 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2155 coloring->F for the coarser grids from the finest 2156 */ 2157 if (coloring->F) { 2158 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2159 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2160 if (m1 != m2) { 2161 coloring->F = 0; 2162 } 2163 } 2164 2165 if (coloring->F) { 2166 w1 = coloring->F; 2167 coloring->F = 0; 2168 } else { 2169 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2170 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2171 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2172 } 2173 2174 /* 2175 Compute all the scale factors and share with other processors 2176 */ 2177 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2178 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2179 for (k=0; k<coloring->ncolors; k++) { 2180 /* 2181 Loop over each column associated with color adding the 2182 perturbation to the vector w3. 2183 */ 2184 for (l=0; l<coloring->ncolumns[k]; l++) { 2185 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2186 dx = xx[col]; 2187 if (dx == 0.0) dx = 1.0; 2188 #if !defined(PETSC_USE_COMPLEX) 2189 if (dx < umin && dx >= 0.0) dx = umin; 2190 else if (dx < 0.0 && dx > -umin) dx = -umin; 2191 #else 2192 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2193 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2194 #endif 2195 dx *= epsilon; 2196 vscale_array[col] = 1.0/dx; 2197 } 2198 } 2199 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2200 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2201 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2202 2203 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2204 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2205 2206 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2207 else vscaleforrow = coloring->columnsforrow; 2208 2209 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2210 /* 2211 Loop over each color 2212 */ 2213 for (k=0; k<coloring->ncolors; k++) { 2214 coloring->currentcolor = k; 2215 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2216 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2217 /* 2218 Loop over each column associated with color adding the 2219 perturbation to the vector w3. 2220 */ 2221 for (l=0; l<coloring->ncolumns[k]; l++) { 2222 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2223 dx = xx[col]; 2224 if (dx == 0.0) dx = 1.0; 2225 #if !defined(PETSC_USE_COMPLEX) 2226 if (dx < umin && dx >= 0.0) dx = umin; 2227 else if (dx < 0.0 && dx > -umin) dx = -umin; 2228 #else 2229 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2230 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2231 #endif 2232 dx *= epsilon; 2233 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2234 w3_array[col] += dx; 2235 } 2236 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2237 2238 /* 2239 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2240 */ 2241 2242 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2243 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2244 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2245 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2246 2247 /* 2248 Loop over rows of vector, putting results into Jacobian matrix 2249 */ 2250 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2251 for (l=0; l<coloring->nrows[k]; l++) { 2252 row = coloring->rows[k][l]; 2253 col = coloring->columnsforrow[k][l]; 2254 y[row] *= vscale_array[vscaleforrow[k][l]]; 2255 srow = row + start; 2256 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2257 } 2258 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2259 } 2260 coloring->currentcolor = k; 2261 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2262 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2263 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2264 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2265 PetscFunctionReturn(0); 2266 } 2267 2268 #undef __FUNCT__ 2269 #define __FUNCT__ "MatAXPY_SeqAIJ" 2270 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2271 { 2272 PetscErrorCode ierr; 2273 PetscInt i; 2274 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2275 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2276 2277 PetscFunctionBegin; 2278 if (str == SAME_NONZERO_PATTERN) { 2279 PetscScalar alpha = a; 2280 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2281 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2282 if (y->xtoy && y->XtoY != X) { 2283 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2284 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2285 } 2286 if (!y->xtoy) { /* get xtoy */ 2287 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2288 y->XtoY = X; 2289 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2290 } 2291 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2292 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2293 } else { 2294 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2295 } 2296 PetscFunctionReturn(0); 2297 } 2298 2299 #undef __FUNCT__ 2300 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2301 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2302 { 2303 PetscErrorCode ierr; 2304 2305 PetscFunctionBegin; 2306 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2307 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2308 PetscFunctionReturn(0); 2309 } 2310 2311 #undef __FUNCT__ 2312 #define __FUNCT__ "MatConjugate_SeqAIJ" 2313 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2314 { 2315 #if defined(PETSC_USE_COMPLEX) 2316 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2317 PetscInt i,nz; 2318 PetscScalar *a; 2319 2320 PetscFunctionBegin; 2321 nz = aij->nz; 2322 a = aij->a; 2323 for (i=0; i<nz; i++) { 2324 a[i] = PetscConj(a[i]); 2325 } 2326 #else 2327 PetscFunctionBegin; 2328 #endif 2329 PetscFunctionReturn(0); 2330 } 2331 2332 #undef __FUNCT__ 2333 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2334 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2335 { 2336 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2337 PetscErrorCode ierr; 2338 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2339 PetscReal atmp; 2340 PetscScalar *x; 2341 MatScalar *aa; 2342 2343 PetscFunctionBegin; 2344 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2345 aa = a->a; 2346 ai = a->i; 2347 aj = a->j; 2348 2349 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2350 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2351 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2352 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2353 for (i=0; i<m; i++) { 2354 ncols = ai[1] - ai[0]; ai++; 2355 x[i] = 0.0; 2356 for (j=0; j<ncols; j++){ 2357 atmp = PetscAbsScalar(*aa); 2358 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2359 aa++; aj++; 2360 } 2361 } 2362 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2363 PetscFunctionReturn(0); 2364 } 2365 2366 #undef __FUNCT__ 2367 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2368 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2369 { 2370 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2371 PetscErrorCode ierr; 2372 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2373 PetscScalar *x; 2374 MatScalar *aa; 2375 2376 PetscFunctionBegin; 2377 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2378 aa = a->a; 2379 ai = a->i; 2380 aj = a->j; 2381 2382 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2383 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2384 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2385 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2386 for (i=0; i<m; i++) { 2387 ncols = ai[1] - ai[0]; ai++; 2388 if (ncols == A->cmap->n) { /* row is dense */ 2389 x[i] = *aa; if (idx) idx[i] = 0; 2390 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2391 x[i] = 0.0; 2392 if (idx) { 2393 idx[i] = 0; /* in case ncols is zero */ 2394 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2395 if (aj[j] > j) { 2396 idx[i] = j; 2397 break; 2398 } 2399 } 2400 } 2401 } 2402 for (j=0; j<ncols; j++){ 2403 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2404 aa++; aj++; 2405 } 2406 } 2407 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2408 PetscFunctionReturn(0); 2409 } 2410 2411 #undef __FUNCT__ 2412 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2413 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2414 { 2415 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2416 PetscErrorCode ierr; 2417 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2418 PetscReal atmp; 2419 PetscScalar *x; 2420 MatScalar *aa; 2421 2422 PetscFunctionBegin; 2423 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2424 aa = a->a; 2425 ai = a->i; 2426 aj = a->j; 2427 2428 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2429 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2430 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2431 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2432 for (i=0; i<m; i++) { 2433 ncols = ai[1] - ai[0]; ai++; 2434 if (ncols) { 2435 /* Get first nonzero */ 2436 for(j = 0; j < ncols; j++) { 2437 atmp = PetscAbsScalar(aa[j]); 2438 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2439 } 2440 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2441 } else { 2442 x[i] = 0.0; if (idx) idx[i] = 0; 2443 } 2444 for(j = 0; j < ncols; j++) { 2445 atmp = PetscAbsScalar(*aa); 2446 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2447 aa++; aj++; 2448 } 2449 } 2450 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2451 PetscFunctionReturn(0); 2452 } 2453 2454 #undef __FUNCT__ 2455 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2456 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2457 { 2458 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2459 PetscErrorCode ierr; 2460 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2461 PetscScalar *x; 2462 MatScalar *aa; 2463 2464 PetscFunctionBegin; 2465 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2466 aa = a->a; 2467 ai = a->i; 2468 aj = a->j; 2469 2470 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2471 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2472 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2473 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2474 for (i=0; i<m; i++) { 2475 ncols = ai[1] - ai[0]; ai++; 2476 if (ncols == A->cmap->n) { /* row is dense */ 2477 x[i] = *aa; if (idx) idx[i] = 0; 2478 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2479 x[i] = 0.0; 2480 if (idx) { /* find first implicit 0.0 in the row */ 2481 idx[i] = 0; /* in case ncols is zero */ 2482 for (j=0;j<ncols;j++) { 2483 if (aj[j] > j) { 2484 idx[i] = j; 2485 break; 2486 } 2487 } 2488 } 2489 } 2490 for (j=0; j<ncols; j++){ 2491 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2492 aa++; aj++; 2493 } 2494 } 2495 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2496 PetscFunctionReturn(0); 2497 } 2498 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2499 /* -------------------------------------------------------------------*/ 2500 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2501 MatGetRow_SeqAIJ, 2502 MatRestoreRow_SeqAIJ, 2503 MatMult_SeqAIJ, 2504 /* 4*/ MatMultAdd_SeqAIJ, 2505 MatMultTranspose_SeqAIJ, 2506 MatMultTransposeAdd_SeqAIJ, 2507 0, 2508 0, 2509 0, 2510 /*10*/ 0, 2511 MatLUFactor_SeqAIJ, 2512 0, 2513 MatSOR_SeqAIJ, 2514 MatTranspose_SeqAIJ, 2515 /*15*/ MatGetInfo_SeqAIJ, 2516 MatEqual_SeqAIJ, 2517 MatGetDiagonal_SeqAIJ, 2518 MatDiagonalScale_SeqAIJ, 2519 MatNorm_SeqAIJ, 2520 /*20*/ 0, 2521 MatAssemblyEnd_SeqAIJ, 2522 MatSetOption_SeqAIJ, 2523 MatZeroEntries_SeqAIJ, 2524 /*24*/ MatZeroRows_SeqAIJ, 2525 0, 2526 0, 2527 0, 2528 0, 2529 /*29*/ MatSetUpPreallocation_SeqAIJ, 2530 0, 2531 0, 2532 MatGetArray_SeqAIJ, 2533 MatRestoreArray_SeqAIJ, 2534 /*34*/ MatDuplicate_SeqAIJ, 2535 0, 2536 0, 2537 MatILUFactor_SeqAIJ, 2538 0, 2539 /*39*/ MatAXPY_SeqAIJ, 2540 MatGetSubMatrices_SeqAIJ, 2541 MatIncreaseOverlap_SeqAIJ, 2542 MatGetValues_SeqAIJ, 2543 MatCopy_SeqAIJ, 2544 /*44*/ MatGetRowMax_SeqAIJ, 2545 MatScale_SeqAIJ, 2546 0, 2547 MatDiagonalSet_SeqAIJ, 2548 0, 2549 /*49*/ MatSetBlockSize_SeqAIJ, 2550 MatGetRowIJ_SeqAIJ, 2551 MatRestoreRowIJ_SeqAIJ, 2552 MatGetColumnIJ_SeqAIJ, 2553 MatRestoreColumnIJ_SeqAIJ, 2554 /*54*/ MatFDColoringCreate_SeqAIJ, 2555 0, 2556 0, 2557 MatPermute_SeqAIJ, 2558 0, 2559 /*59*/ 0, 2560 MatDestroy_SeqAIJ, 2561 MatView_SeqAIJ, 2562 0, 2563 0, 2564 /*64*/ 0, 2565 0, 2566 0, 2567 0, 2568 0, 2569 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2570 MatGetRowMinAbs_SeqAIJ, 2571 0, 2572 MatSetColoring_SeqAIJ, 2573 #if defined(PETSC_HAVE_ADIC) 2574 MatSetValuesAdic_SeqAIJ, 2575 #else 2576 0, 2577 #endif 2578 /*74*/ MatSetValuesAdifor_SeqAIJ, 2579 MatFDColoringApply_AIJ, 2580 0, 2581 0, 2582 0, 2583 /*79*/ 0, 2584 0, 2585 0, 2586 0, 2587 MatLoad_SeqAIJ, 2588 /*84*/ MatIsSymmetric_SeqAIJ, 2589 MatIsHermitian_SeqAIJ, 2590 0, 2591 0, 2592 0, 2593 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2594 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2595 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2596 MatPtAP_Basic, 2597 MatPtAPSymbolic_SeqAIJ, 2598 /*94*/ MatPtAPNumeric_SeqAIJ, 2599 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2600 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2601 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2602 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2603 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2604 0, 2605 0, 2606 MatConjugate_SeqAIJ, 2607 0, 2608 /*104*/MatSetValuesRow_SeqAIJ, 2609 MatRealPart_SeqAIJ, 2610 MatImaginaryPart_SeqAIJ, 2611 0, 2612 0, 2613 /*109*/0, 2614 0, 2615 MatGetRowMin_SeqAIJ, 2616 0, 2617 MatMissingDiagonal_SeqAIJ, 2618 /*114*/0, 2619 0, 2620 0, 2621 0, 2622 0, 2623 /*119*/0, 2624 0, 2625 0, 2626 0 2627 }; 2628 2629 EXTERN_C_BEGIN 2630 #undef __FUNCT__ 2631 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2632 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2633 { 2634 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2635 PetscInt i,nz,n; 2636 2637 PetscFunctionBegin; 2638 2639 nz = aij->maxnz; 2640 n = mat->rmap->n; 2641 for (i=0; i<nz; i++) { 2642 aij->j[i] = indices[i]; 2643 } 2644 aij->nz = nz; 2645 for (i=0; i<n; i++) { 2646 aij->ilen[i] = aij->imax[i]; 2647 } 2648 2649 PetscFunctionReturn(0); 2650 } 2651 EXTERN_C_END 2652 2653 #undef __FUNCT__ 2654 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2655 /*@ 2656 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2657 in the matrix. 2658 2659 Input Parameters: 2660 + mat - the SeqAIJ matrix 2661 - indices - the column indices 2662 2663 Level: advanced 2664 2665 Notes: 2666 This can be called if you have precomputed the nonzero structure of the 2667 matrix and want to provide it to the matrix object to improve the performance 2668 of the MatSetValues() operation. 2669 2670 You MUST have set the correct numbers of nonzeros per row in the call to 2671 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2672 2673 MUST be called before any calls to MatSetValues(); 2674 2675 The indices should start with zero, not one. 2676 2677 @*/ 2678 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2679 { 2680 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2681 2682 PetscFunctionBegin; 2683 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2684 PetscValidPointer(indices,2); 2685 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2686 if (f) { 2687 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2688 } else { 2689 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2690 } 2691 PetscFunctionReturn(0); 2692 } 2693 2694 /* ----------------------------------------------------------------------------------------*/ 2695 2696 EXTERN_C_BEGIN 2697 #undef __FUNCT__ 2698 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2699 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2700 { 2701 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2702 PetscErrorCode ierr; 2703 size_t nz = aij->i[mat->rmap->n]; 2704 2705 PetscFunctionBegin; 2706 if (aij->nonew != 1) { 2707 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2708 } 2709 2710 /* allocate space for values if not already there */ 2711 if (!aij->saved_values) { 2712 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2713 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2714 } 2715 2716 /* copy values over */ 2717 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2718 PetscFunctionReturn(0); 2719 } 2720 EXTERN_C_END 2721 2722 #undef __FUNCT__ 2723 #define __FUNCT__ "MatStoreValues" 2724 /*@ 2725 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2726 example, reuse of the linear part of a Jacobian, while recomputing the 2727 nonlinear portion. 2728 2729 Collect on Mat 2730 2731 Input Parameters: 2732 . mat - the matrix (currently only AIJ matrices support this option) 2733 2734 Level: advanced 2735 2736 Common Usage, with SNESSolve(): 2737 $ Create Jacobian matrix 2738 $ Set linear terms into matrix 2739 $ Apply boundary conditions to matrix, at this time matrix must have 2740 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2741 $ boundary conditions again will not change the nonzero structure 2742 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2743 $ ierr = MatStoreValues(mat); 2744 $ Call SNESSetJacobian() with matrix 2745 $ In your Jacobian routine 2746 $ ierr = MatRetrieveValues(mat); 2747 $ Set nonlinear terms in matrix 2748 2749 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2750 $ // build linear portion of Jacobian 2751 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2752 $ ierr = MatStoreValues(mat); 2753 $ loop over nonlinear iterations 2754 $ ierr = MatRetrieveValues(mat); 2755 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2756 $ // call MatAssemblyBegin/End() on matrix 2757 $ Solve linear system with Jacobian 2758 $ endloop 2759 2760 Notes: 2761 Matrix must already be assemblied before calling this routine 2762 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2763 calling this routine. 2764 2765 When this is called multiple times it overwrites the previous set of stored values 2766 and does not allocated additional space. 2767 2768 .seealso: MatRetrieveValues() 2769 2770 @*/ 2771 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2772 { 2773 PetscErrorCode ierr,(*f)(Mat); 2774 2775 PetscFunctionBegin; 2776 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2777 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2778 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2779 2780 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2781 if (f) { 2782 ierr = (*f)(mat);CHKERRQ(ierr); 2783 } else { 2784 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2785 } 2786 PetscFunctionReturn(0); 2787 } 2788 2789 EXTERN_C_BEGIN 2790 #undef __FUNCT__ 2791 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2792 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2793 { 2794 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2795 PetscErrorCode ierr; 2796 PetscInt nz = aij->i[mat->rmap->n]; 2797 2798 PetscFunctionBegin; 2799 if (aij->nonew != 1) { 2800 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2801 } 2802 if (!aij->saved_values) { 2803 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2804 } 2805 /* copy values over */ 2806 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2807 PetscFunctionReturn(0); 2808 } 2809 EXTERN_C_END 2810 2811 #undef __FUNCT__ 2812 #define __FUNCT__ "MatRetrieveValues" 2813 /*@ 2814 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2815 example, reuse of the linear part of a Jacobian, while recomputing the 2816 nonlinear portion. 2817 2818 Collect on Mat 2819 2820 Input Parameters: 2821 . mat - the matrix (currently on AIJ matrices support this option) 2822 2823 Level: advanced 2824 2825 .seealso: MatStoreValues() 2826 2827 @*/ 2828 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2829 { 2830 PetscErrorCode ierr,(*f)(Mat); 2831 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2834 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2835 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2836 2837 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2838 if (f) { 2839 ierr = (*f)(mat);CHKERRQ(ierr); 2840 } else { 2841 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2842 } 2843 PetscFunctionReturn(0); 2844 } 2845 2846 2847 /* --------------------------------------------------------------------------------*/ 2848 #undef __FUNCT__ 2849 #define __FUNCT__ "MatCreateSeqAIJ" 2850 /*@C 2851 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2852 (the default parallel PETSc format). For good matrix assembly performance 2853 the user should preallocate the matrix storage by setting the parameter nz 2854 (or the array nnz). By setting these parameters accurately, performance 2855 during matrix assembly can be increased by more than a factor of 50. 2856 2857 Collective on MPI_Comm 2858 2859 Input Parameters: 2860 + comm - MPI communicator, set to PETSC_COMM_SELF 2861 . m - number of rows 2862 . n - number of columns 2863 . nz - number of nonzeros per row (same for all rows) 2864 - nnz - array containing the number of nonzeros in the various rows 2865 (possibly different for each row) or PETSC_NULL 2866 2867 Output Parameter: 2868 . A - the matrix 2869 2870 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2871 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2872 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2873 2874 Notes: 2875 If nnz is given then nz is ignored 2876 2877 The AIJ format (also called the Yale sparse matrix format or 2878 compressed row storage), is fully compatible with standard Fortran 77 2879 storage. That is, the stored row and column indices can begin at 2880 either one (as in Fortran) or zero. See the users' manual for details. 2881 2882 Specify the preallocated storage with either nz or nnz (not both). 2883 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2884 allocation. For large problems you MUST preallocate memory or you 2885 will get TERRIBLE performance, see the users' manual chapter on matrices. 2886 2887 By default, this format uses inodes (identical nodes) when possible, to 2888 improve numerical efficiency of matrix-vector products and solves. We 2889 search for consecutive rows with the same nonzero structure, thereby 2890 reusing matrix information to achieve increased efficiency. 2891 2892 Options Database Keys: 2893 + -mat_no_inode - Do not use inodes 2894 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2895 - -mat_aij_oneindex - Internally use indexing starting at 1 2896 rather than 0. Note that when calling MatSetValues(), 2897 the user still MUST index entries starting at 0! 2898 2899 Level: intermediate 2900 2901 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2902 2903 @*/ 2904 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2905 { 2906 PetscErrorCode ierr; 2907 2908 PetscFunctionBegin; 2909 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2910 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2911 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2912 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2913 PetscFunctionReturn(0); 2914 } 2915 2916 #undef __FUNCT__ 2917 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2918 /*@C 2919 MatSeqAIJSetPreallocation - For good matrix assembly performance 2920 the user should preallocate the matrix storage by setting the parameter nz 2921 (or the array nnz). By setting these parameters accurately, performance 2922 during matrix assembly can be increased by more than a factor of 50. 2923 2924 Collective on MPI_Comm 2925 2926 Input Parameters: 2927 + B - The matrix-free 2928 . nz - number of nonzeros per row (same for all rows) 2929 - nnz - array containing the number of nonzeros in the various rows 2930 (possibly different for each row) or PETSC_NULL 2931 2932 Notes: 2933 If nnz is given then nz is ignored 2934 2935 The AIJ format (also called the Yale sparse matrix format or 2936 compressed row storage), is fully compatible with standard Fortran 77 2937 storage. That is, the stored row and column indices can begin at 2938 either one (as in Fortran) or zero. See the users' manual for details. 2939 2940 Specify the preallocated storage with either nz or nnz (not both). 2941 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2942 allocation. For large problems you MUST preallocate memory or you 2943 will get TERRIBLE performance, see the users' manual chapter on matrices. 2944 2945 You can call MatGetInfo() to get information on how effective the preallocation was; 2946 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2947 You can also run with the option -info and look for messages with the string 2948 malloc in them to see if additional memory allocation was needed. 2949 2950 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2951 entries or columns indices 2952 2953 By default, this format uses inodes (identical nodes) when possible, to 2954 improve numerical efficiency of matrix-vector products and solves. We 2955 search for consecutive rows with the same nonzero structure, thereby 2956 reusing matrix information to achieve increased efficiency. 2957 2958 Options Database Keys: 2959 + -mat_no_inode - Do not use inodes 2960 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2961 - -mat_aij_oneindex - Internally use indexing starting at 1 2962 rather than 0. Note that when calling MatSetValues(), 2963 the user still MUST index entries starting at 0! 2964 2965 Level: intermediate 2966 2967 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2968 2969 @*/ 2970 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2971 { 2972 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2973 2974 PetscFunctionBegin; 2975 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2976 if (f) { 2977 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2978 } 2979 PetscFunctionReturn(0); 2980 } 2981 2982 EXTERN_C_BEGIN 2983 #undef __FUNCT__ 2984 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2985 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2986 { 2987 Mat_SeqAIJ *b; 2988 PetscTruth skipallocation = PETSC_FALSE; 2989 PetscErrorCode ierr; 2990 PetscInt i; 2991 2992 PetscFunctionBegin; 2993 2994 if (nz == MAT_SKIP_ALLOCATION) { 2995 skipallocation = PETSC_TRUE; 2996 nz = 0; 2997 } 2998 2999 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3000 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3001 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3002 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3003 3004 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3005 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3006 if (nnz) { 3007 for (i=0; i<B->rmap->n; i++) { 3008 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 3009 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3010 } 3011 } 3012 3013 B->preallocated = PETSC_TRUE; 3014 b = (Mat_SeqAIJ*)B->data; 3015 3016 if (!skipallocation) { 3017 if (!b->imax) { 3018 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3019 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3020 } 3021 if (!nnz) { 3022 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3023 else if (nz <= 0) nz = 1; 3024 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3025 nz = nz*B->rmap->n; 3026 } else { 3027 nz = 0; 3028 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3029 } 3030 /* b->ilen will count nonzeros in each row so far. */ 3031 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3032 3033 /* allocate the matrix space */ 3034 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3035 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3036 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3037 b->i[0] = 0; 3038 for (i=1; i<B->rmap->n+1; i++) { 3039 b->i[i] = b->i[i-1] + b->imax[i-1]; 3040 } 3041 b->singlemalloc = PETSC_TRUE; 3042 b->free_a = PETSC_TRUE; 3043 b->free_ij = PETSC_TRUE; 3044 } else { 3045 b->free_a = PETSC_FALSE; 3046 b->free_ij = PETSC_FALSE; 3047 } 3048 3049 b->nz = 0; 3050 b->maxnz = nz; 3051 B->info.nz_unneeded = (double)b->maxnz; 3052 PetscFunctionReturn(0); 3053 } 3054 EXTERN_C_END 3055 3056 #undef __FUNCT__ 3057 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3058 /*@ 3059 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3060 3061 Input Parameters: 3062 + B - the matrix 3063 . i - the indices into j for the start of each row (starts with zero) 3064 . j - the column indices for each row (starts with zero) these must be sorted for each row 3065 - v - optional values in the matrix 3066 3067 Level: developer 3068 3069 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3070 3071 .keywords: matrix, aij, compressed row, sparse, sequential 3072 3073 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3074 @*/ 3075 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3076 { 3077 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3078 PetscErrorCode ierr; 3079 3080 PetscFunctionBegin; 3081 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3082 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3083 if (f) { 3084 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3085 } 3086 PetscFunctionReturn(0); 3087 } 3088 3089 EXTERN_C_BEGIN 3090 #undef __FUNCT__ 3091 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3092 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3093 { 3094 PetscInt i; 3095 PetscInt m,n; 3096 PetscInt nz; 3097 PetscInt *nnz, nz_max = 0; 3098 PetscScalar *values; 3099 PetscErrorCode ierr; 3100 3101 PetscFunctionBegin; 3102 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3103 3104 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3105 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3106 for(i = 0; i < m; i++) { 3107 nz = Ii[i+1]- Ii[i]; 3108 nz_max = PetscMax(nz_max, nz); 3109 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3110 nnz[i] = nz; 3111 } 3112 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3113 ierr = PetscFree(nnz);CHKERRQ(ierr); 3114 3115 if (v) { 3116 values = (PetscScalar*) v; 3117 } else { 3118 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3119 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3120 } 3121 3122 for(i = 0; i < m; i++) { 3123 nz = Ii[i+1] - Ii[i]; 3124 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3125 } 3126 3127 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3128 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3129 3130 if (!v) { 3131 ierr = PetscFree(values);CHKERRQ(ierr); 3132 } 3133 PetscFunctionReturn(0); 3134 } 3135 EXTERN_C_END 3136 3137 #include "../src/mat/impls/dense/seq/dense.h" 3138 #include "private/petscaxpy.h" 3139 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3142 /* 3143 Computes (B'*A')' since computing B*A directly is untenable 3144 3145 n p p 3146 ( ) ( ) ( ) 3147 m ( A ) * n ( B ) = m ( C ) 3148 ( ) ( ) ( ) 3149 3150 */ 3151 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3152 { 3153 PetscErrorCode ierr; 3154 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3155 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3156 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3157 PetscInt i,n,m,q,p; 3158 const PetscInt *ii,*idx; 3159 const PetscScalar *b,*a,*a_q; 3160 PetscScalar *c,*c_q; 3161 3162 PetscFunctionBegin; 3163 m = A->rmap->n; 3164 n = A->cmap->n; 3165 p = B->cmap->n; 3166 a = sub_a->v; 3167 b = sub_b->a; 3168 c = sub_c->v; 3169 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3170 3171 ii = sub_b->i; 3172 idx = sub_b->j; 3173 for (i=0; i<n; i++) { 3174 q = ii[i+1] - ii[i]; 3175 while (q-->0) { 3176 c_q = c + m*(*idx); 3177 a_q = a + m*i; 3178 PetscAXPY(c_q,*b,a_q,m); 3179 idx++; 3180 b++; 3181 } 3182 } 3183 PetscFunctionReturn(0); 3184 } 3185 3186 #undef __FUNCT__ 3187 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3188 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3189 { 3190 PetscErrorCode ierr; 3191 PetscInt m=A->rmap->n,n=B->cmap->n; 3192 Mat Cmat; 3193 3194 PetscFunctionBegin; 3195 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3196 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3197 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3198 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3199 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3200 Cmat->assembled = PETSC_TRUE; 3201 *C = Cmat; 3202 PetscFunctionReturn(0); 3203 } 3204 3205 /* ----------------------------------------------------------------*/ 3206 #undef __FUNCT__ 3207 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3208 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3209 { 3210 PetscErrorCode ierr; 3211 3212 PetscFunctionBegin; 3213 if (scall == MAT_INITIAL_MATRIX){ 3214 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3215 } 3216 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3217 PetscFunctionReturn(0); 3218 } 3219 3220 3221 /*MC 3222 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3223 based on compressed sparse row format. 3224 3225 Options Database Keys: 3226 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3227 3228 Level: beginner 3229 3230 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3231 M*/ 3232 3233 EXTERN_C_BEGIN 3234 #if defined(PETSC_HAVE_PASTIX) 3235 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3236 #endif 3237 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3238 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3239 #endif 3240 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3241 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3242 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3243 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3244 #if defined(PETSC_HAVE_MUMPS) 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3246 #endif 3247 #if defined(PETSC_HAVE_SUPERLU) 3248 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3249 #endif 3250 #if defined(PETSC_HAVE_SUPERLU_DIST) 3251 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3252 #endif 3253 #if defined(PETSC_HAVE_SPOOLES) 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3255 #endif 3256 #if defined(PETSC_HAVE_UMFPACK) 3257 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3258 #endif 3259 #if defined(PETSC_HAVE_CHOLMOD) 3260 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3261 #endif 3262 #if defined(PETSC_HAVE_LUSOL) 3263 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3264 #endif 3265 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3266 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3267 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3268 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3269 #endif 3270 EXTERN_C_END 3271 3272 3273 EXTERN_C_BEGIN 3274 #undef __FUNCT__ 3275 #define __FUNCT__ "MatCreate_SeqAIJ" 3276 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3277 { 3278 Mat_SeqAIJ *b; 3279 PetscErrorCode ierr; 3280 PetscMPIInt size; 3281 3282 PetscFunctionBegin; 3283 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3284 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3285 3286 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3287 B->data = (void*)b; 3288 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3289 B->mapping = 0; 3290 b->row = 0; 3291 b->col = 0; 3292 b->icol = 0; 3293 b->reallocs = 0; 3294 b->ignorezeroentries = PETSC_FALSE; 3295 b->roworiented = PETSC_TRUE; 3296 b->nonew = 0; 3297 b->diag = 0; 3298 b->solve_work = 0; 3299 B->spptr = 0; 3300 b->saved_values = 0; 3301 b->idiag = 0; 3302 b->mdiag = 0; 3303 b->ssor_work = 0; 3304 b->omega = 1.0; 3305 b->fshift = 0.0; 3306 b->idiagvalid = PETSC_FALSE; 3307 b->keepnonzeropattern = PETSC_FALSE; 3308 b->xtoy = 0; 3309 b->XtoY = 0; 3310 b->compressedrow.use = PETSC_FALSE; 3311 b->compressedrow.nrows = B->rmap->n; 3312 b->compressedrow.i = PETSC_NULL; 3313 b->compressedrow.rindex = PETSC_NULL; 3314 b->compressedrow.checked = PETSC_FALSE; 3315 B->same_nonzero = PETSC_FALSE; 3316 3317 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3318 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3319 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3320 "MatGetFactor_seqaij_matlab", 3321 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3322 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3324 #endif 3325 #if defined(PETSC_HAVE_PASTIX) 3326 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3327 "MatGetFactor_seqaij_pastix", 3328 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3329 #endif 3330 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3332 "MatGetFactor_seqaij_essl", 3333 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3334 #endif 3335 #if defined(PETSC_HAVE_SUPERLU) 3336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3337 "MatGetFactor_seqaij_superlu", 3338 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3339 #endif 3340 #if defined(PETSC_HAVE_SUPERLU_DIST) 3341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3342 "MatGetFactor_seqaij_superlu_dist", 3343 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3344 #endif 3345 #if defined(PETSC_HAVE_SPOOLES) 3346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3347 "MatGetFactor_seqaij_spooles", 3348 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3349 #endif 3350 #if defined(PETSC_HAVE_MUMPS) 3351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3352 "MatGetFactor_aij_mumps", 3353 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3354 #endif 3355 #if defined(PETSC_HAVE_UMFPACK) 3356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3357 "MatGetFactor_seqaij_umfpack", 3358 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3359 #endif 3360 #if defined(PETSC_HAVE_CHOLMOD) 3361 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3362 "MatGetFactor_seqaij_cholmod", 3363 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3364 #endif 3365 #if defined(PETSC_HAVE_LUSOL) 3366 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3367 "MatGetFactor_seqaij_lusol", 3368 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3369 #endif 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3371 "MatGetFactor_seqaij_petsc", 3372 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3373 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3374 "MatGetFactorAvailable_seqaij_petsc", 3375 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3376 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3377 "MatGetFactor_seqaij_bas", 3378 MatGetFactor_seqaij_bas); 3379 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3380 "MatSeqAIJSetColumnIndices_SeqAIJ", 3381 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3383 "MatStoreValues_SeqAIJ", 3384 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3386 "MatRetrieveValues_SeqAIJ", 3387 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3388 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3389 "MatConvert_SeqAIJ_SeqSBAIJ", 3390 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3391 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3392 "MatConvert_SeqAIJ_SeqBAIJ", 3393 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3394 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3395 "MatConvert_SeqAIJ_SeqCSRPERM", 3396 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3397 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3398 "MatConvert_SeqAIJ_SeqCRL", 3399 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3400 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3401 "MatIsTranspose_SeqAIJ", 3402 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3404 "MatIsHermitianTranspose_SeqAIJ", 3405 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3406 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3407 "MatSeqAIJSetPreallocation_SeqAIJ", 3408 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3410 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3411 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3412 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3413 "MatReorderForNonzeroDiagonal_SeqAIJ", 3414 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3415 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3416 "MatMatMult_SeqDense_SeqAIJ", 3417 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3418 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3419 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3420 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3421 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3422 "MatMatMultNumeric_SeqDense_SeqAIJ", 3423 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3424 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3425 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3426 PetscFunctionReturn(0); 3427 } 3428 EXTERN_C_END 3429 3430 #undef __FUNCT__ 3431 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3432 /* 3433 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3434 */ 3435 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3436 { 3437 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3438 PetscErrorCode ierr; 3439 PetscInt i,m = A->rmap->n; 3440 3441 PetscFunctionBegin; 3442 c = (Mat_SeqAIJ*)C->data; 3443 3444 C->factortype = A->factortype; 3445 c->row = 0; 3446 c->col = 0; 3447 c->icol = 0; 3448 c->reallocs = 0; 3449 3450 C->assembled = PETSC_TRUE; 3451 3452 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3453 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3454 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3455 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3456 3457 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3458 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3459 for (i=0; i<m; i++) { 3460 c->imax[i] = a->imax[i]; 3461 c->ilen[i] = a->ilen[i]; 3462 } 3463 3464 /* allocate the matrix space */ 3465 if (mallocmatspace){ 3466 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3467 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3468 c->singlemalloc = PETSC_TRUE; 3469 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3470 if (m > 0) { 3471 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3472 if (cpvalues == MAT_COPY_VALUES) { 3473 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3474 } else { 3475 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3476 } 3477 } 3478 } 3479 3480 c->ignorezeroentries = a->ignorezeroentries; 3481 c->roworiented = a->roworiented; 3482 c->nonew = a->nonew; 3483 if (a->diag) { 3484 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3485 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3486 for (i=0; i<m; i++) { 3487 c->diag[i] = a->diag[i]; 3488 } 3489 } else c->diag = 0; 3490 c->solve_work = 0; 3491 c->saved_values = 0; 3492 c->idiag = 0; 3493 c->ssor_work = 0; 3494 c->keepnonzeropattern = a->keepnonzeropattern; 3495 c->free_a = PETSC_TRUE; 3496 c->free_ij = PETSC_TRUE; 3497 c->xtoy = 0; 3498 c->XtoY = 0; 3499 3500 c->nz = a->nz; 3501 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3502 C->preallocated = PETSC_TRUE; 3503 3504 c->compressedrow.use = a->compressedrow.use; 3505 c->compressedrow.nrows = a->compressedrow.nrows; 3506 c->compressedrow.checked = a->compressedrow.checked; 3507 if (a->compressedrow.checked && a->compressedrow.use){ 3508 i = a->compressedrow.nrows; 3509 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3510 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3511 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3512 } else { 3513 c->compressedrow.use = PETSC_FALSE; 3514 c->compressedrow.i = PETSC_NULL; 3515 c->compressedrow.rindex = PETSC_NULL; 3516 } 3517 C->same_nonzero = A->same_nonzero; 3518 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3519 3520 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3521 PetscFunctionReturn(0); 3522 } 3523 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3526 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3527 { 3528 PetscErrorCode ierr; 3529 3530 PetscFunctionBegin; 3531 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3532 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3533 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3534 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3535 PetscFunctionReturn(0); 3536 } 3537 3538 #undef __FUNCT__ 3539 #define __FUNCT__ "MatLoad_SeqAIJ" 3540 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,Mat newMat) 3541 { 3542 Mat_SeqAIJ *a; 3543 PetscErrorCode ierr; 3544 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3545 int fd; 3546 PetscMPIInt size; 3547 MPI_Comm comm; 3548 3549 PetscFunctionBegin; 3550 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3551 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3552 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3553 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3554 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3555 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3556 M = header[1]; N = header[2]; nz = header[3]; 3557 3558 if (nz < 0) { 3559 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3560 } 3561 3562 /* read in row lengths */ 3563 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3564 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3565 3566 /* check if sum of rowlengths is same as nz */ 3567 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3568 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3569 3570 /* set global size if not set already*/ 3571 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3572 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3573 } else { 3574 /* if sizes and type are already set, check if the vector global sizes are correct */ 3575 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3576 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 3577 } 3578 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3579 a = (Mat_SeqAIJ*)newMat->data; 3580 3581 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3582 3583 /* read in nonzero values */ 3584 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3585 3586 /* set matrix "i" values */ 3587 a->i[0] = 0; 3588 for (i=1; i<= M; i++) { 3589 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3590 a->ilen[i-1] = rowlengths[i-1]; 3591 } 3592 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3593 3594 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3595 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3596 PetscFunctionReturn(0); 3597 } 3598 3599 #undef __FUNCT__ 3600 #define __FUNCT__ "MatEqual_SeqAIJ" 3601 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3602 { 3603 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3604 PetscErrorCode ierr; 3605 #if defined(PETSC_USE_COMPLEX) 3606 PetscInt k; 3607 #endif 3608 3609 PetscFunctionBegin; 3610 /* If the matrix dimensions are not equal,or no of nonzeros */ 3611 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3612 *flg = PETSC_FALSE; 3613 PetscFunctionReturn(0); 3614 } 3615 3616 /* if the a->i are the same */ 3617 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3618 if (!*flg) PetscFunctionReturn(0); 3619 3620 /* if a->j are the same */ 3621 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3622 if (!*flg) PetscFunctionReturn(0); 3623 3624 /* if a->a are the same */ 3625 #if defined(PETSC_USE_COMPLEX) 3626 for (k=0; k<a->nz; k++){ 3627 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3628 *flg = PETSC_FALSE; 3629 PetscFunctionReturn(0); 3630 } 3631 } 3632 #else 3633 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3634 #endif 3635 PetscFunctionReturn(0); 3636 } 3637 3638 #undef __FUNCT__ 3639 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3640 /*@ 3641 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3642 provided by the user. 3643 3644 Collective on MPI_Comm 3645 3646 Input Parameters: 3647 + comm - must be an MPI communicator of size 1 3648 . m - number of rows 3649 . n - number of columns 3650 . i - row indices 3651 . j - column indices 3652 - a - matrix values 3653 3654 Output Parameter: 3655 . mat - the matrix 3656 3657 Level: intermediate 3658 3659 Notes: 3660 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3661 once the matrix is destroyed 3662 3663 You cannot set new nonzero locations into this matrix, that will generate an error. 3664 3665 The i and j indices are 0 based 3666 3667 The format which is used for the sparse matrix input, is equivalent to a 3668 row-major ordering.. i.e for the following matrix, the input data expected is 3669 as shown: 3670 3671 1 0 0 3672 2 0 3 3673 4 5 6 3674 3675 i = {0,1,3,6} [size = nrow+1 = 3+1] 3676 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3677 v = {1,2,3,4,5,6} [size = nz = 6] 3678 3679 3680 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3681 3682 @*/ 3683 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3684 { 3685 PetscErrorCode ierr; 3686 PetscInt ii; 3687 Mat_SeqAIJ *aij; 3688 #if defined(PETSC_USE_DEBUG) 3689 PetscInt jj; 3690 #endif 3691 3692 PetscFunctionBegin; 3693 if (i[0]) { 3694 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3695 } 3696 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3697 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3698 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3699 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3700 aij = (Mat_SeqAIJ*)(*mat)->data; 3701 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3702 3703 aij->i = i; 3704 aij->j = j; 3705 aij->a = a; 3706 aij->singlemalloc = PETSC_FALSE; 3707 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3708 aij->free_a = PETSC_FALSE; 3709 aij->free_ij = PETSC_FALSE; 3710 3711 for (ii=0; ii<m; ii++) { 3712 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3713 #if defined(PETSC_USE_DEBUG) 3714 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3715 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3716 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 3717 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 3718 } 3719 #endif 3720 } 3721 #if defined(PETSC_USE_DEBUG) 3722 for (ii=0; ii<aij->i[m]; ii++) { 3723 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3724 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3725 } 3726 #endif 3727 3728 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3729 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3730 PetscFunctionReturn(0); 3731 } 3732 3733 #undef __FUNCT__ 3734 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3735 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3736 { 3737 PetscErrorCode ierr; 3738 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3739 3740 PetscFunctionBegin; 3741 if (coloring->ctype == IS_COLORING_GLOBAL) { 3742 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3743 a->coloring = coloring; 3744 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3745 PetscInt i,*larray; 3746 ISColoring ocoloring; 3747 ISColoringValue *colors; 3748 3749 /* set coloring for diagonal portion */ 3750 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3751 for (i=0; i<A->cmap->n; i++) { 3752 larray[i] = i; 3753 } 3754 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3755 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3756 for (i=0; i<A->cmap->n; i++) { 3757 colors[i] = coloring->colors[larray[i]]; 3758 } 3759 ierr = PetscFree(larray);CHKERRQ(ierr); 3760 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3761 a->coloring = ocoloring; 3762 } 3763 PetscFunctionReturn(0); 3764 } 3765 3766 #if defined(PETSC_HAVE_ADIC) 3767 EXTERN_C_BEGIN 3768 #include "adic/ad_utils.h" 3769 EXTERN_C_END 3770 3771 #undef __FUNCT__ 3772 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3773 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3774 { 3775 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3776 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3777 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3778 ISColoringValue *color; 3779 3780 PetscFunctionBegin; 3781 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3782 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3783 color = a->coloring->colors; 3784 /* loop over rows */ 3785 for (i=0; i<m; i++) { 3786 nz = ii[i+1] - ii[i]; 3787 /* loop over columns putting computed value into matrix */ 3788 for (j=0; j<nz; j++) { 3789 *v++ = values[color[*jj++]]; 3790 } 3791 values += nlen; /* jump to next row of derivatives */ 3792 } 3793 PetscFunctionReturn(0); 3794 } 3795 #endif 3796 3797 #undef __FUNCT__ 3798 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3799 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3800 { 3801 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3802 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3803 MatScalar *v = a->a; 3804 PetscScalar *values = (PetscScalar *)advalues; 3805 ISColoringValue *color; 3806 3807 PetscFunctionBegin; 3808 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3809 color = a->coloring->colors; 3810 /* loop over rows */ 3811 for (i=0; i<m; i++) { 3812 nz = ii[i+1] - ii[i]; 3813 /* loop over columns putting computed value into matrix */ 3814 for (j=0; j<nz; j++) { 3815 *v++ = values[color[*jj++]]; 3816 } 3817 values += nl; /* jump to next row of derivatives */ 3818 } 3819 PetscFunctionReturn(0); 3820 } 3821 3822 /* 3823 Special version for direct calls from Fortran 3824 */ 3825 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3826 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3827 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3828 #define matsetvaluesseqaij_ matsetvaluesseqaij 3829 #endif 3830 3831 /* Change these macros so can be used in void function */ 3832 #undef CHKERRQ 3833 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3834 #undef SETERRQ2 3835 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3836 3837 EXTERN_C_BEGIN 3838 #undef __FUNCT__ 3839 #define __FUNCT__ "matsetvaluesseqaij_" 3840 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3841 { 3842 Mat A = *AA; 3843 PetscInt m = *mm, n = *nn; 3844 InsertMode is = *isis; 3845 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3846 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3847 PetscInt *imax,*ai,*ailen; 3848 PetscErrorCode ierr; 3849 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3850 MatScalar *ap,value,*aa; 3851 PetscTruth ignorezeroentries = a->ignorezeroentries; 3852 PetscTruth roworiented = a->roworiented; 3853 3854 PetscFunctionBegin; 3855 ierr = MatPreallocated(A);CHKERRQ(ierr); 3856 imax = a->imax; 3857 ai = a->i; 3858 ailen = a->ilen; 3859 aj = a->j; 3860 aa = a->a; 3861 3862 for (k=0; k<m; k++) { /* loop over added rows */ 3863 row = im[k]; 3864 if (row < 0) continue; 3865 #if defined(PETSC_USE_DEBUG) 3866 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3867 #endif 3868 rp = aj + ai[row]; ap = aa + ai[row]; 3869 rmax = imax[row]; nrow = ailen[row]; 3870 low = 0; 3871 high = nrow; 3872 for (l=0; l<n; l++) { /* loop over added columns */ 3873 if (in[l] < 0) continue; 3874 #if defined(PETSC_USE_DEBUG) 3875 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3876 #endif 3877 col = in[l]; 3878 if (roworiented) { 3879 value = v[l + k*n]; 3880 } else { 3881 value = v[k + l*m]; 3882 } 3883 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3884 3885 if (col <= lastcol) low = 0; else high = nrow; 3886 lastcol = col; 3887 while (high-low > 5) { 3888 t = (low+high)/2; 3889 if (rp[t] > col) high = t; 3890 else low = t; 3891 } 3892 for (i=low; i<high; i++) { 3893 if (rp[i] > col) break; 3894 if (rp[i] == col) { 3895 if (is == ADD_VALUES) ap[i] += value; 3896 else ap[i] = value; 3897 goto noinsert; 3898 } 3899 } 3900 if (value == 0.0 && ignorezeroentries) goto noinsert; 3901 if (nonew == 1) goto noinsert; 3902 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3903 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3904 N = nrow++ - 1; a->nz++; high++; 3905 /* shift up all the later entries in this row */ 3906 for (ii=N; ii>=i; ii--) { 3907 rp[ii+1] = rp[ii]; 3908 ap[ii+1] = ap[ii]; 3909 } 3910 rp[i] = col; 3911 ap[i] = value; 3912 noinsert:; 3913 low = i + 1; 3914 } 3915 ailen[row] = nrow; 3916 } 3917 A->same_nonzero = PETSC_FALSE; 3918 PetscFunctionReturnVoid(); 3919 } 3920 EXTERN_C_END 3921