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 2896 Level: intermediate 2897 2898 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2899 2900 @*/ 2901 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2902 { 2903 PetscErrorCode ierr; 2904 2905 PetscFunctionBegin; 2906 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2907 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2908 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2909 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2910 PetscFunctionReturn(0); 2911 } 2912 2913 #undef __FUNCT__ 2914 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2915 /*@C 2916 MatSeqAIJSetPreallocation - For good matrix assembly performance 2917 the user should preallocate the matrix storage by setting the parameter nz 2918 (or the array nnz). By setting these parameters accurately, performance 2919 during matrix assembly can be increased by more than a factor of 50. 2920 2921 Collective on MPI_Comm 2922 2923 Input Parameters: 2924 + B - The matrix-free 2925 . nz - number of nonzeros per row (same for all rows) 2926 - nnz - array containing the number of nonzeros in the various rows 2927 (possibly different for each row) or PETSC_NULL 2928 2929 Notes: 2930 If nnz is given then nz is ignored 2931 2932 The AIJ format (also called the Yale sparse matrix format or 2933 compressed row storage), is fully compatible with standard Fortran 77 2934 storage. That is, the stored row and column indices can begin at 2935 either one (as in Fortran) or zero. See the users' manual for details. 2936 2937 Specify the preallocated storage with either nz or nnz (not both). 2938 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2939 allocation. For large problems you MUST preallocate memory or you 2940 will get TERRIBLE performance, see the users' manual chapter on matrices. 2941 2942 You can call MatGetInfo() to get information on how effective the preallocation was; 2943 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2944 You can also run with the option -info and look for messages with the string 2945 malloc in them to see if additional memory allocation was needed. 2946 2947 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2948 entries or columns indices 2949 2950 By default, this format uses inodes (identical nodes) when possible, to 2951 improve numerical efficiency of matrix-vector products and solves. We 2952 search for consecutive rows with the same nonzero structure, thereby 2953 reusing matrix information to achieve increased efficiency. 2954 2955 Options Database Keys: 2956 + -mat_no_inode - Do not use inodes 2957 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2958 - -mat_aij_oneindex - Internally use indexing starting at 1 2959 rather than 0. Note that when calling MatSetValues(), 2960 the user still MUST index entries starting at 0! 2961 2962 Level: intermediate 2963 2964 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2965 2966 @*/ 2967 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2968 { 2969 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2970 2971 PetscFunctionBegin; 2972 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2973 if (f) { 2974 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2975 } 2976 PetscFunctionReturn(0); 2977 } 2978 2979 EXTERN_C_BEGIN 2980 #undef __FUNCT__ 2981 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2982 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2983 { 2984 Mat_SeqAIJ *b; 2985 PetscTruth skipallocation = PETSC_FALSE; 2986 PetscErrorCode ierr; 2987 PetscInt i; 2988 2989 PetscFunctionBegin; 2990 2991 if (nz == MAT_SKIP_ALLOCATION) { 2992 skipallocation = PETSC_TRUE; 2993 nz = 0; 2994 } 2995 2996 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2997 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2998 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2999 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3000 3001 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3002 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3003 if (nnz) { 3004 for (i=0; i<B->rmap->n; i++) { 3005 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]); 3006 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); 3007 } 3008 } 3009 3010 B->preallocated = PETSC_TRUE; 3011 b = (Mat_SeqAIJ*)B->data; 3012 3013 if (!skipallocation) { 3014 if (!b->imax) { 3015 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3016 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3017 } 3018 if (!nnz) { 3019 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3020 else if (nz <= 0) nz = 1; 3021 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3022 nz = nz*B->rmap->n; 3023 } else { 3024 nz = 0; 3025 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3026 } 3027 /* b->ilen will count nonzeros in each row so far. */ 3028 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3029 3030 /* allocate the matrix space */ 3031 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3032 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3033 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3034 b->i[0] = 0; 3035 for (i=1; i<B->rmap->n+1; i++) { 3036 b->i[i] = b->i[i-1] + b->imax[i-1]; 3037 } 3038 b->singlemalloc = PETSC_TRUE; 3039 b->free_a = PETSC_TRUE; 3040 b->free_ij = PETSC_TRUE; 3041 } else { 3042 b->free_a = PETSC_FALSE; 3043 b->free_ij = PETSC_FALSE; 3044 } 3045 3046 b->nz = 0; 3047 b->maxnz = nz; 3048 B->info.nz_unneeded = (double)b->maxnz; 3049 PetscFunctionReturn(0); 3050 } 3051 EXTERN_C_END 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3055 /*@ 3056 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3057 3058 Input Parameters: 3059 + B - the matrix 3060 . i - the indices into j for the start of each row (starts with zero) 3061 . j - the column indices for each row (starts with zero) these must be sorted for each row 3062 - v - optional values in the matrix 3063 3064 Level: developer 3065 3066 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3067 3068 .keywords: matrix, aij, compressed row, sparse, sequential 3069 3070 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3071 @*/ 3072 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3073 { 3074 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3075 PetscErrorCode ierr; 3076 3077 PetscFunctionBegin; 3078 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3079 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3080 if (f) { 3081 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3082 } 3083 PetscFunctionReturn(0); 3084 } 3085 3086 EXTERN_C_BEGIN 3087 #undef __FUNCT__ 3088 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3089 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3090 { 3091 PetscInt i; 3092 PetscInt m,n; 3093 PetscInt nz; 3094 PetscInt *nnz, nz_max = 0; 3095 PetscScalar *values; 3096 PetscErrorCode ierr; 3097 3098 PetscFunctionBegin; 3099 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3100 3101 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3102 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3103 for(i = 0; i < m; i++) { 3104 nz = Ii[i+1]- Ii[i]; 3105 nz_max = PetscMax(nz_max, nz); 3106 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3107 nnz[i] = nz; 3108 } 3109 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3110 ierr = PetscFree(nnz);CHKERRQ(ierr); 3111 3112 if (v) { 3113 values = (PetscScalar*) v; 3114 } else { 3115 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3116 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3117 } 3118 3119 for(i = 0; i < m; i++) { 3120 nz = Ii[i+1] - Ii[i]; 3121 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3122 } 3123 3124 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3125 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3126 3127 if (!v) { 3128 ierr = PetscFree(values);CHKERRQ(ierr); 3129 } 3130 PetscFunctionReturn(0); 3131 } 3132 EXTERN_C_END 3133 3134 #include "../src/mat/impls/dense/seq/dense.h" 3135 #include "private/petscaxpy.h" 3136 3137 #undef __FUNCT__ 3138 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3139 /* 3140 Computes (B'*A')' since computing B*A directly is untenable 3141 3142 n p p 3143 ( ) ( ) ( ) 3144 m ( A ) * n ( B ) = m ( C ) 3145 ( ) ( ) ( ) 3146 3147 */ 3148 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3149 { 3150 PetscErrorCode ierr; 3151 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3152 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3153 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3154 PetscInt i,n,m,q,p; 3155 const PetscInt *ii,*idx; 3156 const PetscScalar *b,*a,*a_q; 3157 PetscScalar *c,*c_q; 3158 3159 PetscFunctionBegin; 3160 m = A->rmap->n; 3161 n = A->cmap->n; 3162 p = B->cmap->n; 3163 a = sub_a->v; 3164 b = sub_b->a; 3165 c = sub_c->v; 3166 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3167 3168 ii = sub_b->i; 3169 idx = sub_b->j; 3170 for (i=0; i<n; i++) { 3171 q = ii[i+1] - ii[i]; 3172 while (q-->0) { 3173 c_q = c + m*(*idx); 3174 a_q = a + m*i; 3175 PetscAXPY(c_q,*b,a_q,m); 3176 idx++; 3177 b++; 3178 } 3179 } 3180 PetscFunctionReturn(0); 3181 } 3182 3183 #undef __FUNCT__ 3184 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3185 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3186 { 3187 PetscErrorCode ierr; 3188 PetscInt m=A->rmap->n,n=B->cmap->n; 3189 Mat Cmat; 3190 3191 PetscFunctionBegin; 3192 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); 3193 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3194 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3195 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3196 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3197 Cmat->assembled = PETSC_TRUE; 3198 *C = Cmat; 3199 PetscFunctionReturn(0); 3200 } 3201 3202 /* ----------------------------------------------------------------*/ 3203 #undef __FUNCT__ 3204 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3205 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3206 { 3207 PetscErrorCode ierr; 3208 3209 PetscFunctionBegin; 3210 if (scall == MAT_INITIAL_MATRIX){ 3211 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3212 } 3213 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3214 PetscFunctionReturn(0); 3215 } 3216 3217 3218 /*MC 3219 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3220 based on compressed sparse row format. 3221 3222 Options Database Keys: 3223 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3224 3225 Level: beginner 3226 3227 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3228 M*/ 3229 3230 EXTERN_C_BEGIN 3231 #if defined(PETSC_HAVE_PASTIX) 3232 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3233 #endif 3234 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3235 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3236 #endif 3237 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3238 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3239 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3240 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3241 #if defined(PETSC_HAVE_MUMPS) 3242 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3243 #endif 3244 #if defined(PETSC_HAVE_SUPERLU) 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3246 #endif 3247 #if defined(PETSC_HAVE_SUPERLU_DIST) 3248 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3249 #endif 3250 #if defined(PETSC_HAVE_SPOOLES) 3251 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3252 #endif 3253 #if defined(PETSC_HAVE_UMFPACK) 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3255 #endif 3256 #if defined(PETSC_HAVE_CHOLMOD) 3257 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3258 #endif 3259 #if defined(PETSC_HAVE_LUSOL) 3260 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3261 #endif 3262 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3263 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3264 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3265 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3266 #endif 3267 EXTERN_C_END 3268 3269 3270 EXTERN_C_BEGIN 3271 #undef __FUNCT__ 3272 #define __FUNCT__ "MatCreate_SeqAIJ" 3273 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3274 { 3275 Mat_SeqAIJ *b; 3276 PetscErrorCode ierr; 3277 PetscMPIInt size; 3278 3279 PetscFunctionBegin; 3280 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3281 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3282 3283 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3284 B->data = (void*)b; 3285 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3286 B->mapping = 0; 3287 b->row = 0; 3288 b->col = 0; 3289 b->icol = 0; 3290 b->reallocs = 0; 3291 b->ignorezeroentries = PETSC_FALSE; 3292 b->roworiented = PETSC_TRUE; 3293 b->nonew = 0; 3294 b->diag = 0; 3295 b->solve_work = 0; 3296 B->spptr = 0; 3297 b->saved_values = 0; 3298 b->idiag = 0; 3299 b->mdiag = 0; 3300 b->ssor_work = 0; 3301 b->omega = 1.0; 3302 b->fshift = 0.0; 3303 b->idiagvalid = PETSC_FALSE; 3304 b->keepnonzeropattern = PETSC_FALSE; 3305 b->xtoy = 0; 3306 b->XtoY = 0; 3307 b->compressedrow.use = PETSC_FALSE; 3308 b->compressedrow.nrows = B->rmap->n; 3309 b->compressedrow.i = PETSC_NULL; 3310 b->compressedrow.rindex = PETSC_NULL; 3311 b->compressedrow.checked = PETSC_FALSE; 3312 B->same_nonzero = PETSC_FALSE; 3313 3314 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3315 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3316 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3317 "MatGetFactor_seqaij_matlab", 3318 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3319 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3320 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3321 #endif 3322 #if defined(PETSC_HAVE_PASTIX) 3323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3324 "MatGetFactor_seqaij_pastix", 3325 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3326 #endif 3327 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3329 "MatGetFactor_seqaij_essl", 3330 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3331 #endif 3332 #if defined(PETSC_HAVE_SUPERLU) 3333 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3334 "MatGetFactor_seqaij_superlu", 3335 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3336 #endif 3337 #if defined(PETSC_HAVE_SUPERLU_DIST) 3338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3339 "MatGetFactor_seqaij_superlu_dist", 3340 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3341 #endif 3342 #if defined(PETSC_HAVE_SPOOLES) 3343 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3344 "MatGetFactor_seqaij_spooles", 3345 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3346 #endif 3347 #if defined(PETSC_HAVE_MUMPS) 3348 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3349 "MatGetFactor_aij_mumps", 3350 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3351 #endif 3352 #if defined(PETSC_HAVE_UMFPACK) 3353 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3354 "MatGetFactor_seqaij_umfpack", 3355 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3356 #endif 3357 #if defined(PETSC_HAVE_CHOLMOD) 3358 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3359 "MatGetFactor_seqaij_cholmod", 3360 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3361 #endif 3362 #if defined(PETSC_HAVE_LUSOL) 3363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3364 "MatGetFactor_seqaij_lusol", 3365 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3366 #endif 3367 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3368 "MatGetFactor_seqaij_petsc", 3369 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3371 "MatGetFactorAvailable_seqaij_petsc", 3372 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3373 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3374 "MatGetFactor_seqaij_bas", 3375 MatGetFactor_seqaij_bas); 3376 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3377 "MatSeqAIJSetColumnIndices_SeqAIJ", 3378 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3379 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3380 "MatStoreValues_SeqAIJ", 3381 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3383 "MatRetrieveValues_SeqAIJ", 3384 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3386 "MatConvert_SeqAIJ_SeqSBAIJ", 3387 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3388 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3389 "MatConvert_SeqAIJ_SeqBAIJ", 3390 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3391 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3392 "MatConvert_SeqAIJ_SeqCSRPERM", 3393 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3394 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3395 "MatConvert_SeqAIJ_SeqCRL", 3396 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3397 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3398 "MatIsTranspose_SeqAIJ", 3399 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3400 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3401 "MatIsHermitianTranspose_SeqAIJ", 3402 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3404 "MatSeqAIJSetPreallocation_SeqAIJ", 3405 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3406 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3407 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3408 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3410 "MatReorderForNonzeroDiagonal_SeqAIJ", 3411 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3412 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3413 "MatMatMult_SeqDense_SeqAIJ", 3414 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3415 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3416 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3417 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3418 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3419 "MatMatMultNumeric_SeqDense_SeqAIJ", 3420 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3421 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3422 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3423 PetscFunctionReturn(0); 3424 } 3425 EXTERN_C_END 3426 3427 #undef __FUNCT__ 3428 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3429 /* 3430 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3431 */ 3432 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3433 { 3434 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3435 PetscErrorCode ierr; 3436 PetscInt i,m = A->rmap->n; 3437 3438 PetscFunctionBegin; 3439 c = (Mat_SeqAIJ*)C->data; 3440 3441 C->factortype = A->factortype; 3442 c->row = 0; 3443 c->col = 0; 3444 c->icol = 0; 3445 c->reallocs = 0; 3446 3447 C->assembled = PETSC_TRUE; 3448 3449 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3450 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3451 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3452 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3453 3454 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3455 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3456 for (i=0; i<m; i++) { 3457 c->imax[i] = a->imax[i]; 3458 c->ilen[i] = a->ilen[i]; 3459 } 3460 3461 /* allocate the matrix space */ 3462 if (mallocmatspace){ 3463 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3464 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3465 c->singlemalloc = PETSC_TRUE; 3466 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3467 if (m > 0) { 3468 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3469 if (cpvalues == MAT_COPY_VALUES) { 3470 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3471 } else { 3472 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3473 } 3474 } 3475 } 3476 3477 c->ignorezeroentries = a->ignorezeroentries; 3478 c->roworiented = a->roworiented; 3479 c->nonew = a->nonew; 3480 if (a->diag) { 3481 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3482 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3483 for (i=0; i<m; i++) { 3484 c->diag[i] = a->diag[i]; 3485 } 3486 } else c->diag = 0; 3487 c->solve_work = 0; 3488 c->saved_values = 0; 3489 c->idiag = 0; 3490 c->ssor_work = 0; 3491 c->keepnonzeropattern = a->keepnonzeropattern; 3492 c->free_a = PETSC_TRUE; 3493 c->free_ij = PETSC_TRUE; 3494 c->xtoy = 0; 3495 c->XtoY = 0; 3496 3497 c->nz = a->nz; 3498 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3499 C->preallocated = PETSC_TRUE; 3500 3501 c->compressedrow.use = a->compressedrow.use; 3502 c->compressedrow.nrows = a->compressedrow.nrows; 3503 c->compressedrow.checked = a->compressedrow.checked; 3504 if (a->compressedrow.checked && a->compressedrow.use){ 3505 i = a->compressedrow.nrows; 3506 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3507 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3508 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3509 } else { 3510 c->compressedrow.use = PETSC_FALSE; 3511 c->compressedrow.i = PETSC_NULL; 3512 c->compressedrow.rindex = PETSC_NULL; 3513 } 3514 C->same_nonzero = A->same_nonzero; 3515 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3516 3517 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3518 PetscFunctionReturn(0); 3519 } 3520 3521 #undef __FUNCT__ 3522 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3523 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3524 { 3525 PetscErrorCode ierr; 3526 3527 PetscFunctionBegin; 3528 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3529 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3530 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3531 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3532 PetscFunctionReturn(0); 3533 } 3534 3535 #undef __FUNCT__ 3536 #define __FUNCT__ "MatLoad_SeqAIJ" 3537 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3538 { 3539 Mat_SeqAIJ *a; 3540 PetscErrorCode ierr; 3541 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3542 int fd; 3543 PetscMPIInt size; 3544 MPI_Comm comm; 3545 3546 PetscFunctionBegin; 3547 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3548 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3549 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3550 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3551 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3552 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3553 M = header[1]; N = header[2]; nz = header[3]; 3554 3555 if (nz < 0) { 3556 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3557 } 3558 3559 /* read in row lengths */ 3560 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3561 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3562 3563 /* check if sum of rowlengths is same as nz */ 3564 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3565 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); 3566 3567 /* set global size if not set already*/ 3568 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3569 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3570 } else { 3571 /* if sizes and type are already set, check if the vector global sizes are correct */ 3572 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3573 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); 3574 } 3575 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3576 a = (Mat_SeqAIJ*)newMat->data; 3577 3578 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3579 3580 /* read in nonzero values */ 3581 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3582 3583 /* set matrix "i" values */ 3584 a->i[0] = 0; 3585 for (i=1; i<= M; i++) { 3586 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3587 a->ilen[i-1] = rowlengths[i-1]; 3588 } 3589 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3590 3591 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3592 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3593 PetscFunctionReturn(0); 3594 } 3595 3596 #undef __FUNCT__ 3597 #define __FUNCT__ "MatEqual_SeqAIJ" 3598 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3599 { 3600 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3601 PetscErrorCode ierr; 3602 #if defined(PETSC_USE_COMPLEX) 3603 PetscInt k; 3604 #endif 3605 3606 PetscFunctionBegin; 3607 /* If the matrix dimensions are not equal,or no of nonzeros */ 3608 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3609 *flg = PETSC_FALSE; 3610 PetscFunctionReturn(0); 3611 } 3612 3613 /* if the a->i are the same */ 3614 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3615 if (!*flg) PetscFunctionReturn(0); 3616 3617 /* if a->j are the same */ 3618 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3619 if (!*flg) PetscFunctionReturn(0); 3620 3621 /* if a->a are the same */ 3622 #if defined(PETSC_USE_COMPLEX) 3623 for (k=0; k<a->nz; k++){ 3624 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3625 *flg = PETSC_FALSE; 3626 PetscFunctionReturn(0); 3627 } 3628 } 3629 #else 3630 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3631 #endif 3632 PetscFunctionReturn(0); 3633 } 3634 3635 #undef __FUNCT__ 3636 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3637 /*@ 3638 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3639 provided by the user. 3640 3641 Collective on MPI_Comm 3642 3643 Input Parameters: 3644 + comm - must be an MPI communicator of size 1 3645 . m - number of rows 3646 . n - number of columns 3647 . i - row indices 3648 . j - column indices 3649 - a - matrix values 3650 3651 Output Parameter: 3652 . mat - the matrix 3653 3654 Level: intermediate 3655 3656 Notes: 3657 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3658 once the matrix is destroyed 3659 3660 You cannot set new nonzero locations into this matrix, that will generate an error. 3661 3662 The i and j indices are 0 based 3663 3664 The format which is used for the sparse matrix input, is equivalent to a 3665 row-major ordering.. i.e for the following matrix, the input data expected is 3666 as shown: 3667 3668 1 0 0 3669 2 0 3 3670 4 5 6 3671 3672 i = {0,1,3,6} [size = nrow+1 = 3+1] 3673 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3674 v = {1,2,3,4,5,6} [size = nz = 6] 3675 3676 3677 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3678 3679 @*/ 3680 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3681 { 3682 PetscErrorCode ierr; 3683 PetscInt ii; 3684 Mat_SeqAIJ *aij; 3685 #if defined(PETSC_USE_DEBUG) 3686 PetscInt jj; 3687 #endif 3688 3689 PetscFunctionBegin; 3690 if (i[0]) { 3691 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3692 } 3693 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3694 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3695 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3696 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3697 aij = (Mat_SeqAIJ*)(*mat)->data; 3698 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3699 3700 aij->i = i; 3701 aij->j = j; 3702 aij->a = a; 3703 aij->singlemalloc = PETSC_FALSE; 3704 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3705 aij->free_a = PETSC_FALSE; 3706 aij->free_ij = PETSC_FALSE; 3707 3708 for (ii=0; ii<m; ii++) { 3709 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3710 #if defined(PETSC_USE_DEBUG) 3711 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]); 3712 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3713 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); 3714 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); 3715 } 3716 #endif 3717 } 3718 #if defined(PETSC_USE_DEBUG) 3719 for (ii=0; ii<aij->i[m]; ii++) { 3720 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3721 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]); 3722 } 3723 #endif 3724 3725 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3726 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3727 PetscFunctionReturn(0); 3728 } 3729 3730 #undef __FUNCT__ 3731 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3732 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3733 { 3734 PetscErrorCode ierr; 3735 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3736 3737 PetscFunctionBegin; 3738 if (coloring->ctype == IS_COLORING_GLOBAL) { 3739 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3740 a->coloring = coloring; 3741 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3742 PetscInt i,*larray; 3743 ISColoring ocoloring; 3744 ISColoringValue *colors; 3745 3746 /* set coloring for diagonal portion */ 3747 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3748 for (i=0; i<A->cmap->n; i++) { 3749 larray[i] = i; 3750 } 3751 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3752 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3753 for (i=0; i<A->cmap->n; i++) { 3754 colors[i] = coloring->colors[larray[i]]; 3755 } 3756 ierr = PetscFree(larray);CHKERRQ(ierr); 3757 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3758 a->coloring = ocoloring; 3759 } 3760 PetscFunctionReturn(0); 3761 } 3762 3763 #if defined(PETSC_HAVE_ADIC) 3764 EXTERN_C_BEGIN 3765 #include "adic/ad_utils.h" 3766 EXTERN_C_END 3767 3768 #undef __FUNCT__ 3769 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3770 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3771 { 3772 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3773 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3774 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3775 ISColoringValue *color; 3776 3777 PetscFunctionBegin; 3778 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3779 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3780 color = a->coloring->colors; 3781 /* loop over rows */ 3782 for (i=0; i<m; i++) { 3783 nz = ii[i+1] - ii[i]; 3784 /* loop over columns putting computed value into matrix */ 3785 for (j=0; j<nz; j++) { 3786 *v++ = values[color[*jj++]]; 3787 } 3788 values += nlen; /* jump to next row of derivatives */ 3789 } 3790 PetscFunctionReturn(0); 3791 } 3792 #endif 3793 3794 #undef __FUNCT__ 3795 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3796 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3797 { 3798 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3799 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3800 MatScalar *v = a->a; 3801 PetscScalar *values = (PetscScalar *)advalues; 3802 ISColoringValue *color; 3803 3804 PetscFunctionBegin; 3805 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3806 color = a->coloring->colors; 3807 /* loop over rows */ 3808 for (i=0; i<m; i++) { 3809 nz = ii[i+1] - ii[i]; 3810 /* loop over columns putting computed value into matrix */ 3811 for (j=0; j<nz; j++) { 3812 *v++ = values[color[*jj++]]; 3813 } 3814 values += nl; /* jump to next row of derivatives */ 3815 } 3816 PetscFunctionReturn(0); 3817 } 3818 3819 /* 3820 Special version for direct calls from Fortran 3821 */ 3822 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3823 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3824 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3825 #define matsetvaluesseqaij_ matsetvaluesseqaij 3826 #endif 3827 3828 /* Change these macros so can be used in void function */ 3829 #undef CHKERRQ 3830 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3831 #undef SETERRQ2 3832 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3833 3834 EXTERN_C_BEGIN 3835 #undef __FUNCT__ 3836 #define __FUNCT__ "matsetvaluesseqaij_" 3837 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3838 { 3839 Mat A = *AA; 3840 PetscInt m = *mm, n = *nn; 3841 InsertMode is = *isis; 3842 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3843 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3844 PetscInt *imax,*ai,*ailen; 3845 PetscErrorCode ierr; 3846 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3847 MatScalar *ap,value,*aa; 3848 PetscTruth ignorezeroentries = a->ignorezeroentries; 3849 PetscTruth roworiented = a->roworiented; 3850 3851 PetscFunctionBegin; 3852 ierr = MatPreallocated(A);CHKERRQ(ierr); 3853 imax = a->imax; 3854 ai = a->i; 3855 ailen = a->ilen; 3856 aj = a->j; 3857 aa = a->a; 3858 3859 for (k=0; k<m; k++) { /* loop over added rows */ 3860 row = im[k]; 3861 if (row < 0) continue; 3862 #if defined(PETSC_USE_DEBUG) 3863 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3864 #endif 3865 rp = aj + ai[row]; ap = aa + ai[row]; 3866 rmax = imax[row]; nrow = ailen[row]; 3867 low = 0; 3868 high = nrow; 3869 for (l=0; l<n; l++) { /* loop over added columns */ 3870 if (in[l] < 0) continue; 3871 #if defined(PETSC_USE_DEBUG) 3872 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3873 #endif 3874 col = in[l]; 3875 if (roworiented) { 3876 value = v[l + k*n]; 3877 } else { 3878 value = v[k + l*m]; 3879 } 3880 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3881 3882 if (col <= lastcol) low = 0; else high = nrow; 3883 lastcol = col; 3884 while (high-low > 5) { 3885 t = (low+high)/2; 3886 if (rp[t] > col) high = t; 3887 else low = t; 3888 } 3889 for (i=low; i<high; i++) { 3890 if (rp[i] > col) break; 3891 if (rp[i] == col) { 3892 if (is == ADD_VALUES) ap[i] += value; 3893 else ap[i] = value; 3894 goto noinsert; 3895 } 3896 } 3897 if (value == 0.0 && ignorezeroentries) goto noinsert; 3898 if (nonew == 1) goto noinsert; 3899 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3900 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3901 N = nrow++ - 1; a->nz++; high++; 3902 /* shift up all the later entries in this row */ 3903 for (ii=N; ii>=i; ii--) { 3904 rp[ii+1] = rp[ii]; 3905 ap[ii+1] = ap[ii]; 3906 } 3907 rp[i] = col; 3908 ap[i] = value; 3909 noinsert:; 3910 low = i + 1; 3911 } 3912 ailen[row] = nrow; 3913 } 3914 A->same_nonzero = PETSC_FALSE; 3915 PetscFunctionReturnVoid(); 3916 } 3917 EXTERN_C_END 3918