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_seqaijperm_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 #if defined(PETSC_HAVE_CUDA) 1111 /* 1112 ierr = VecView(xx,0);CHKERRQ(ierr); 1113 ierr = VecView(zz,0);CHKERRQ(ierr); 1114 ierr = MatView(A,0);CHKERRQ(ierr); 1115 */ 1116 #endif 1117 PetscFunctionReturn(0); 1118 } 1119 1120 /* 1121 Adds diagonal pointers to sparse matrix structure. 1122 */ 1123 #undef __FUNCT__ 1124 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1125 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1126 { 1127 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1128 PetscErrorCode ierr; 1129 PetscInt i,j,m = A->rmap->n; 1130 1131 PetscFunctionBegin; 1132 if (!a->diag) { 1133 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1134 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1135 } 1136 for (i=0; i<A->rmap->n; i++) { 1137 a->diag[i] = a->i[i+1]; 1138 for (j=a->i[i]; j<a->i[i+1]; j++) { 1139 if (a->j[j] == i) { 1140 a->diag[i] = j; 1141 break; 1142 } 1143 } 1144 } 1145 PetscFunctionReturn(0); 1146 } 1147 1148 /* 1149 Checks for missing diagonals 1150 */ 1151 #undef __FUNCT__ 1152 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1153 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1154 { 1155 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1156 PetscInt *diag,*jj = a->j,i; 1157 1158 PetscFunctionBegin; 1159 *missing = PETSC_FALSE; 1160 if (A->rmap->n > 0 && !jj) { 1161 *missing = PETSC_TRUE; 1162 if (d) *d = 0; 1163 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1164 } else { 1165 diag = a->diag; 1166 for (i=0; i<A->rmap->n; i++) { 1167 if (jj[diag[i]] != i) { 1168 *missing = PETSC_TRUE; 1169 if (d) *d = i; 1170 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1171 } 1172 } 1173 } 1174 PetscFunctionReturn(0); 1175 } 1176 1177 EXTERN_C_BEGIN 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1180 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1181 { 1182 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1183 PetscErrorCode ierr; 1184 PetscInt i,*diag,m = A->rmap->n; 1185 MatScalar *v = a->a; 1186 PetscScalar *idiag,*mdiag; 1187 1188 PetscFunctionBegin; 1189 if (a->idiagvalid) PetscFunctionReturn(0); 1190 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1191 diag = a->diag; 1192 if (!a->idiag) { 1193 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1194 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1195 v = a->a; 1196 } 1197 mdiag = a->mdiag; 1198 idiag = a->idiag; 1199 1200 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1201 for (i=0; i<m; i++) { 1202 mdiag[i] = v[diag[i]]; 1203 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1204 idiag[i] = 1.0/v[diag[i]]; 1205 } 1206 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1207 } else { 1208 for (i=0; i<m; i++) { 1209 mdiag[i] = v[diag[i]]; 1210 idiag[i] = omega/(fshift + v[diag[i]]); 1211 } 1212 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1213 } 1214 a->idiagvalid = PETSC_TRUE; 1215 PetscFunctionReturn(0); 1216 } 1217 EXTERN_C_END 1218 1219 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "MatSOR_SeqAIJ" 1222 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1223 { 1224 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1225 PetscScalar *x,d,sum,*t,scale; 1226 const MatScalar *v = a->a,*idiag=0,*mdiag; 1227 const PetscScalar *b, *bs,*xb, *ts; 1228 PetscErrorCode ierr; 1229 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1230 const PetscInt *idx,*diag; 1231 1232 PetscFunctionBegin; 1233 its = its*lits; 1234 1235 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1236 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1237 a->fshift = fshift; 1238 a->omega = omega; 1239 1240 diag = a->diag; 1241 t = a->ssor_work; 1242 idiag = a->idiag; 1243 mdiag = a->mdiag; 1244 1245 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1246 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1247 CHKMEMQ; 1248 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1249 if (flag == SOR_APPLY_UPPER) { 1250 /* apply (U + D/omega) to the vector */ 1251 bs = b; 1252 for (i=0; i<m; i++) { 1253 d = fshift + mdiag[i]; 1254 n = a->i[i+1] - diag[i] - 1; 1255 idx = a->j + diag[i] + 1; 1256 v = a->a + diag[i] + 1; 1257 sum = b[i]*d/omega; 1258 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1259 x[i] = sum; 1260 } 1261 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1262 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1263 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1264 PetscFunctionReturn(0); 1265 } 1266 1267 if (flag == SOR_APPLY_LOWER) { 1268 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1269 } else if (flag & SOR_EISENSTAT) { 1270 /* Let A = L + U + D; where L is lower trianglar, 1271 U is upper triangular, E = D/omega; This routine applies 1272 1273 (L + E)^{-1} A (U + E)^{-1} 1274 1275 to a vector efficiently using Eisenstat's trick. 1276 */ 1277 scale = (2.0/omega) - 1.0; 1278 1279 /* x = (E + U)^{-1} b */ 1280 for (i=m-1; i>=0; i--) { 1281 n = a->i[i+1] - diag[i] - 1; 1282 idx = a->j + diag[i] + 1; 1283 v = a->a + diag[i] + 1; 1284 sum = b[i]; 1285 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1286 x[i] = sum*idiag[i]; 1287 } 1288 1289 /* t = b - (2*E - D)x */ 1290 v = a->a; 1291 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1292 1293 /* t = (E + L)^{-1}t */ 1294 ts = t; 1295 diag = a->diag; 1296 for (i=0; i<m; i++) { 1297 n = diag[i] - a->i[i]; 1298 idx = a->j + a->i[i]; 1299 v = a->a + a->i[i]; 1300 sum = t[i]; 1301 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1302 t[i] = sum*idiag[i]; 1303 /* x = x + t */ 1304 x[i] += t[i]; 1305 } 1306 1307 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1308 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1309 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1310 PetscFunctionReturn(0); 1311 } 1312 if (flag & SOR_ZERO_INITIAL_GUESS) { 1313 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1314 for (i=0; i<m; i++) { 1315 n = diag[i] - a->i[i]; 1316 idx = a->j + a->i[i]; 1317 v = a->a + a->i[i]; 1318 sum = b[i]; 1319 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1320 t[i] = sum; 1321 x[i] = sum*idiag[i]; 1322 } 1323 xb = t; 1324 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1325 } else xb = b; 1326 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1327 for (i=m-1; i>=0; i--) { 1328 n = a->i[i+1] - diag[i] - 1; 1329 idx = a->j + diag[i] + 1; 1330 v = a->a + diag[i] + 1; 1331 sum = xb[i]; 1332 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1333 if (xb == b) { 1334 x[i] = sum*idiag[i]; 1335 } else { 1336 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1337 } 1338 } 1339 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1340 } 1341 its--; 1342 } 1343 while (its--) { 1344 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1345 for (i=0; i<m; i++) { 1346 n = a->i[i+1] - a->i[i]; 1347 idx = a->j + a->i[i]; 1348 v = a->a + a->i[i]; 1349 sum = b[i]; 1350 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1351 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1352 } 1353 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1354 } 1355 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1356 for (i=m-1; i>=0; i--) { 1357 n = a->i[i+1] - a->i[i]; 1358 idx = a->j + a->i[i]; 1359 v = a->a + a->i[i]; 1360 sum = b[i]; 1361 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1362 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1363 } 1364 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1365 } 1366 } 1367 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1368 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1369 CHKMEMQ; PetscFunctionReturn(0); 1370 } 1371 1372 1373 #undef __FUNCT__ 1374 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1375 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1376 { 1377 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1378 1379 PetscFunctionBegin; 1380 info->block_size = 1.0; 1381 info->nz_allocated = (double)a->maxnz; 1382 info->nz_used = (double)a->nz; 1383 info->nz_unneeded = (double)(a->maxnz - a->nz); 1384 info->assemblies = (double)A->num_ass; 1385 info->mallocs = (double)A->info.mallocs; 1386 info->memory = ((PetscObject)A)->mem; 1387 if (A->factortype) { 1388 info->fill_ratio_given = A->info.fill_ratio_given; 1389 info->fill_ratio_needed = A->info.fill_ratio_needed; 1390 info->factor_mallocs = A->info.factor_mallocs; 1391 } else { 1392 info->fill_ratio_given = 0; 1393 info->fill_ratio_needed = 0; 1394 info->factor_mallocs = 0; 1395 } 1396 PetscFunctionReturn(0); 1397 } 1398 1399 #undef __FUNCT__ 1400 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1401 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1402 { 1403 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1404 PetscInt i,m = A->rmap->n - 1,d = 0; 1405 PetscErrorCode ierr; 1406 PetscTruth missing; 1407 1408 PetscFunctionBegin; 1409 if (a->keepnonzeropattern) { 1410 for (i=0; i<N; i++) { 1411 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1412 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1413 } 1414 if (diag != 0.0) { 1415 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1416 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1417 for (i=0; i<N; i++) { 1418 a->a[a->diag[rows[i]]] = diag; 1419 } 1420 } 1421 A->same_nonzero = PETSC_TRUE; 1422 } else { 1423 if (diag != 0.0) { 1424 for (i=0; i<N; i++) { 1425 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1426 if (a->ilen[rows[i]] > 0) { 1427 a->ilen[rows[i]] = 1; 1428 a->a[a->i[rows[i]]] = diag; 1429 a->j[a->i[rows[i]]] = rows[i]; 1430 } else { /* in case row was completely empty */ 1431 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1432 } 1433 } 1434 } else { 1435 for (i=0; i<N; i++) { 1436 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1437 a->ilen[rows[i]] = 0; 1438 } 1439 } 1440 A->same_nonzero = PETSC_FALSE; 1441 } 1442 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1443 PetscFunctionReturn(0); 1444 } 1445 1446 #undef __FUNCT__ 1447 #define __FUNCT__ "MatGetRow_SeqAIJ" 1448 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1449 { 1450 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1451 PetscInt *itmp; 1452 1453 PetscFunctionBegin; 1454 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1455 1456 *nz = a->i[row+1] - a->i[row]; 1457 if (v) *v = a->a + a->i[row]; 1458 if (idx) { 1459 itmp = a->j + a->i[row]; 1460 if (*nz) { 1461 *idx = itmp; 1462 } 1463 else *idx = 0; 1464 } 1465 PetscFunctionReturn(0); 1466 } 1467 1468 /* remove this function? */ 1469 #undef __FUNCT__ 1470 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1471 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1472 { 1473 PetscFunctionBegin; 1474 PetscFunctionReturn(0); 1475 } 1476 1477 #undef __FUNCT__ 1478 #define __FUNCT__ "MatNorm_SeqAIJ" 1479 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1480 { 1481 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1482 MatScalar *v = a->a; 1483 PetscReal sum = 0.0; 1484 PetscErrorCode ierr; 1485 PetscInt i,j; 1486 1487 PetscFunctionBegin; 1488 if (type == NORM_FROBENIUS) { 1489 for (i=0; i<a->nz; i++) { 1490 #if defined(PETSC_USE_COMPLEX) 1491 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1492 #else 1493 sum += (*v)*(*v); v++; 1494 #endif 1495 } 1496 *nrm = sqrt(sum); 1497 } else if (type == NORM_1) { 1498 PetscReal *tmp; 1499 PetscInt *jj = a->j; 1500 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1501 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1502 *nrm = 0.0; 1503 for (j=0; j<a->nz; j++) { 1504 tmp[*jj++] += PetscAbsScalar(*v); v++; 1505 } 1506 for (j=0; j<A->cmap->n; j++) { 1507 if (tmp[j] > *nrm) *nrm = tmp[j]; 1508 } 1509 ierr = PetscFree(tmp);CHKERRQ(ierr); 1510 } else if (type == NORM_INFINITY) { 1511 *nrm = 0.0; 1512 for (j=0; j<A->rmap->n; j++) { 1513 v = a->a + a->i[j]; 1514 sum = 0.0; 1515 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1516 sum += PetscAbsScalar(*v); v++; 1517 } 1518 if (sum > *nrm) *nrm = sum; 1519 } 1520 } else { 1521 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1522 } 1523 PetscFunctionReturn(0); 1524 } 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "MatTranspose_SeqAIJ" 1528 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1529 { 1530 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1531 Mat C; 1532 PetscErrorCode ierr; 1533 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1534 MatScalar *array = a->a; 1535 1536 PetscFunctionBegin; 1537 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"); 1538 1539 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1540 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1541 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1542 1543 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1544 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1545 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1546 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1547 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1548 ierr = PetscFree(col);CHKERRQ(ierr); 1549 } else { 1550 C = *B; 1551 } 1552 1553 for (i=0; i<m; i++) { 1554 len = ai[i+1]-ai[i]; 1555 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1556 array += len; 1557 aj += len; 1558 } 1559 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1560 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1561 1562 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1563 *B = C; 1564 } else { 1565 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1566 } 1567 PetscFunctionReturn(0); 1568 } 1569 1570 EXTERN_C_BEGIN 1571 #undef __FUNCT__ 1572 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1573 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1574 { 1575 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1576 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1577 MatScalar *va,*vb; 1578 PetscErrorCode ierr; 1579 PetscInt ma,na,mb,nb, i; 1580 1581 PetscFunctionBegin; 1582 bij = (Mat_SeqAIJ *) B->data; 1583 1584 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1585 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1586 if (ma!=nb || na!=mb){ 1587 *f = PETSC_FALSE; 1588 PetscFunctionReturn(0); 1589 } 1590 aii = aij->i; bii = bij->i; 1591 adx = aij->j; bdx = bij->j; 1592 va = aij->a; vb = bij->a; 1593 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1594 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1595 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1596 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1597 1598 *f = PETSC_TRUE; 1599 for (i=0; i<ma; i++) { 1600 while (aptr[i]<aii[i+1]) { 1601 PetscInt idc,idr; 1602 PetscScalar vc,vr; 1603 /* column/row index/value */ 1604 idc = adx[aptr[i]]; 1605 idr = bdx[bptr[idc]]; 1606 vc = va[aptr[i]]; 1607 vr = vb[bptr[idc]]; 1608 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1609 *f = PETSC_FALSE; 1610 goto done; 1611 } else { 1612 aptr[i]++; 1613 if (B || i!=idc) bptr[idc]++; 1614 } 1615 } 1616 } 1617 done: 1618 ierr = PetscFree(aptr);CHKERRQ(ierr); 1619 if (B) { 1620 ierr = PetscFree(bptr);CHKERRQ(ierr); 1621 } 1622 PetscFunctionReturn(0); 1623 } 1624 EXTERN_C_END 1625 1626 EXTERN_C_BEGIN 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1629 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1630 { 1631 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1632 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1633 MatScalar *va,*vb; 1634 PetscErrorCode ierr; 1635 PetscInt ma,na,mb,nb, i; 1636 1637 PetscFunctionBegin; 1638 bij = (Mat_SeqAIJ *) B->data; 1639 1640 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1641 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1642 if (ma!=nb || na!=mb){ 1643 *f = PETSC_FALSE; 1644 PetscFunctionReturn(0); 1645 } 1646 aii = aij->i; bii = bij->i; 1647 adx = aij->j; bdx = bij->j; 1648 va = aij->a; vb = bij->a; 1649 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1650 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1651 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1652 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1653 1654 *f = PETSC_TRUE; 1655 for (i=0; i<ma; i++) { 1656 while (aptr[i]<aii[i+1]) { 1657 PetscInt idc,idr; 1658 PetscScalar vc,vr; 1659 /* column/row index/value */ 1660 idc = adx[aptr[i]]; 1661 idr = bdx[bptr[idc]]; 1662 vc = va[aptr[i]]; 1663 vr = vb[bptr[idc]]; 1664 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1665 *f = PETSC_FALSE; 1666 goto done; 1667 } else { 1668 aptr[i]++; 1669 if (B || i!=idc) bptr[idc]++; 1670 } 1671 } 1672 } 1673 done: 1674 ierr = PetscFree(aptr);CHKERRQ(ierr); 1675 if (B) { 1676 ierr = PetscFree(bptr);CHKERRQ(ierr); 1677 } 1678 PetscFunctionReturn(0); 1679 } 1680 EXTERN_C_END 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1684 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1685 { 1686 PetscErrorCode ierr; 1687 PetscFunctionBegin; 1688 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 #undef __FUNCT__ 1693 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1694 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1695 { 1696 PetscErrorCode ierr; 1697 PetscFunctionBegin; 1698 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1699 PetscFunctionReturn(0); 1700 } 1701 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1704 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1705 { 1706 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1707 PetscScalar *l,*r,x; 1708 MatScalar *v; 1709 PetscErrorCode ierr; 1710 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1711 1712 PetscFunctionBegin; 1713 if (ll) { 1714 /* The local size is used so that VecMPI can be passed to this routine 1715 by MatDiagonalScale_MPIAIJ */ 1716 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1717 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1718 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1719 v = a->a; 1720 for (i=0; i<m; i++) { 1721 x = l[i]; 1722 M = a->i[i+1] - a->i[i]; 1723 for (j=0; j<M; j++) { (*v++) *= x;} 1724 } 1725 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1726 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1727 } 1728 if (rr) { 1729 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1730 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1731 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1732 v = a->a; jj = a->j; 1733 for (i=0; i<nz; i++) { 1734 (*v++) *= r[*jj++]; 1735 } 1736 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1737 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1738 } 1739 PetscFunctionReturn(0); 1740 } 1741 1742 #undef __FUNCT__ 1743 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1744 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1745 { 1746 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1747 PetscErrorCode ierr; 1748 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1749 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1750 const PetscInt *irow,*icol; 1751 PetscInt nrows,ncols; 1752 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1753 MatScalar *a_new,*mat_a; 1754 Mat C; 1755 PetscTruth stride,sorted; 1756 1757 PetscFunctionBegin; 1758 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1759 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1760 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1761 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1762 1763 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1764 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1765 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1766 1767 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1768 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1769 if (stride && step == 1) { 1770 /* special case of contiguous rows */ 1771 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1772 /* loop over new rows determining lens and starting points */ 1773 for (i=0; i<nrows; i++) { 1774 kstart = ai[irow[i]]; 1775 kend = kstart + ailen[irow[i]]; 1776 for (k=kstart; k<kend; k++) { 1777 if (aj[k] >= first) { 1778 starts[i] = k; 1779 break; 1780 } 1781 } 1782 sum = 0; 1783 while (k < kend) { 1784 if (aj[k++] >= first+ncols) break; 1785 sum++; 1786 } 1787 lens[i] = sum; 1788 } 1789 /* create submatrix */ 1790 if (scall == MAT_REUSE_MATRIX) { 1791 PetscInt n_cols,n_rows; 1792 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1793 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1794 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1795 C = *B; 1796 } else { 1797 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1798 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1799 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1800 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1801 } 1802 c = (Mat_SeqAIJ*)C->data; 1803 1804 /* loop over rows inserting into submatrix */ 1805 a_new = c->a; 1806 j_new = c->j; 1807 i_new = c->i; 1808 1809 for (i=0; i<nrows; i++) { 1810 ii = starts[i]; 1811 lensi = lens[i]; 1812 for (k=0; k<lensi; k++) { 1813 *j_new++ = aj[ii+k] - first; 1814 } 1815 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1816 a_new += lensi; 1817 i_new[i+1] = i_new[i] + lensi; 1818 c->ilen[i] = lensi; 1819 } 1820 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1821 } else { 1822 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1823 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1824 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1825 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1826 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1827 /* determine lens of each row */ 1828 for (i=0; i<nrows; i++) { 1829 kstart = ai[irow[i]]; 1830 kend = kstart + a->ilen[irow[i]]; 1831 lens[i] = 0; 1832 for (k=kstart; k<kend; k++) { 1833 if (smap[aj[k]]) { 1834 lens[i]++; 1835 } 1836 } 1837 } 1838 /* Create and fill new matrix */ 1839 if (scall == MAT_REUSE_MATRIX) { 1840 PetscTruth equal; 1841 1842 c = (Mat_SeqAIJ *)((*B)->data); 1843 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1844 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1845 if (!equal) { 1846 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1847 } 1848 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1849 C = *B; 1850 } else { 1851 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1852 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1853 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1854 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1855 } 1856 c = (Mat_SeqAIJ *)(C->data); 1857 for (i=0; i<nrows; i++) { 1858 row = irow[i]; 1859 kstart = ai[row]; 1860 kend = kstart + a->ilen[row]; 1861 mat_i = c->i[i]; 1862 mat_j = c->j + mat_i; 1863 mat_a = c->a + mat_i; 1864 mat_ilen = c->ilen + i; 1865 for (k=kstart; k<kend; k++) { 1866 if ((tcol=smap[a->j[k]])) { 1867 *mat_j++ = tcol - 1; 1868 *mat_a++ = a->a[k]; 1869 (*mat_ilen)++; 1870 1871 } 1872 } 1873 } 1874 /* Free work space */ 1875 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1876 ierr = PetscFree(smap);CHKERRQ(ierr); 1877 ierr = PetscFree(lens);CHKERRQ(ierr); 1878 } 1879 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1880 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1881 1882 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1883 *B = C; 1884 PetscFunctionReturn(0); 1885 } 1886 1887 #undef __FUNCT__ 1888 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1889 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1890 { 1891 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1892 PetscErrorCode ierr; 1893 Mat outA; 1894 PetscTruth row_identity,col_identity; 1895 1896 PetscFunctionBegin; 1897 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1898 1899 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1900 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1901 1902 outA = inA; 1903 outA->factortype = MAT_FACTOR_LU; 1904 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1905 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1906 a->row = row; 1907 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1908 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1909 a->col = col; 1910 1911 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1912 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1913 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1914 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1915 1916 if (!a->solve_work) { /* this matrix may have been factored before */ 1917 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1918 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1919 } 1920 1921 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1922 if (row_identity && col_identity) { 1923 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1924 } else { 1925 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1926 } 1927 PetscFunctionReturn(0); 1928 } 1929 1930 #undef __FUNCT__ 1931 #define __FUNCT__ "MatScale_SeqAIJ" 1932 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1933 { 1934 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1935 PetscScalar oalpha = alpha; 1936 PetscErrorCode ierr; 1937 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1938 1939 PetscFunctionBegin; 1940 BLASscal_(&bnz,&oalpha,a->a,&one); 1941 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 #undef __FUNCT__ 1946 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1947 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1948 { 1949 PetscErrorCode ierr; 1950 PetscInt i; 1951 1952 PetscFunctionBegin; 1953 if (scall == MAT_INITIAL_MATRIX) { 1954 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1955 } 1956 1957 for (i=0; i<n; i++) { 1958 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1959 } 1960 PetscFunctionReturn(0); 1961 } 1962 1963 #undef __FUNCT__ 1964 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1965 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1966 { 1967 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1968 PetscErrorCode ierr; 1969 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1970 const PetscInt *idx; 1971 PetscInt start,end,*ai,*aj; 1972 PetscBT table; 1973 1974 PetscFunctionBegin; 1975 m = A->rmap->n; 1976 ai = a->i; 1977 aj = a->j; 1978 1979 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1980 1981 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1982 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1983 1984 for (i=0; i<is_max; i++) { 1985 /* Initialize the two local arrays */ 1986 isz = 0; 1987 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1988 1989 /* Extract the indices, assume there can be duplicate entries */ 1990 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1991 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1992 1993 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1994 for (j=0; j<n ; ++j){ 1995 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1996 } 1997 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1998 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1999 2000 k = 0; 2001 for (j=0; j<ov; j++){ /* for each overlap */ 2002 n = isz; 2003 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2004 row = nidx[k]; 2005 start = ai[row]; 2006 end = ai[row+1]; 2007 for (l = start; l<end ; l++){ 2008 val = aj[l] ; 2009 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2010 } 2011 } 2012 } 2013 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2014 } 2015 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2016 ierr = PetscFree(nidx);CHKERRQ(ierr); 2017 PetscFunctionReturn(0); 2018 } 2019 2020 /* -------------------------------------------------------------- */ 2021 #undef __FUNCT__ 2022 #define __FUNCT__ "MatPermute_SeqAIJ" 2023 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2024 { 2025 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2026 PetscErrorCode ierr; 2027 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2028 const PetscInt *row,*col; 2029 PetscInt *cnew,j,*lens; 2030 IS icolp,irowp; 2031 PetscInt *cwork = PETSC_NULL; 2032 PetscScalar *vwork = PETSC_NULL; 2033 2034 PetscFunctionBegin; 2035 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2036 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2037 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2038 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2039 2040 /* determine lengths of permuted rows */ 2041 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2042 for (i=0; i<m; i++) { 2043 lens[row[i]] = a->i[i+1] - a->i[i]; 2044 } 2045 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2046 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2047 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2048 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2049 ierr = PetscFree(lens);CHKERRQ(ierr); 2050 2051 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2052 for (i=0; i<m; i++) { 2053 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2054 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2055 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2056 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2057 } 2058 ierr = PetscFree(cnew);CHKERRQ(ierr); 2059 (*B)->assembled = PETSC_FALSE; 2060 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2061 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2062 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2063 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2064 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2065 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #undef __FUNCT__ 2070 #define __FUNCT__ "MatCopy_SeqAIJ" 2071 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2072 { 2073 PetscErrorCode ierr; 2074 2075 PetscFunctionBegin; 2076 /* If the two matrices have the same copy implementation, use fast copy. */ 2077 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2078 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2079 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2080 2081 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2082 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2083 } 2084 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2085 } else { 2086 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2087 } 2088 PetscFunctionReturn(0); 2089 } 2090 2091 #undef __FUNCT__ 2092 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2093 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2094 { 2095 PetscErrorCode ierr; 2096 2097 PetscFunctionBegin; 2098 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2099 PetscFunctionReturn(0); 2100 } 2101 2102 #undef __FUNCT__ 2103 #define __FUNCT__ "MatGetArray_SeqAIJ" 2104 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2105 { 2106 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2107 PetscFunctionBegin; 2108 *array = a->a; 2109 PetscFunctionReturn(0); 2110 } 2111 2112 #undef __FUNCT__ 2113 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2114 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2115 { 2116 PetscFunctionBegin; 2117 PetscFunctionReturn(0); 2118 } 2119 2120 #undef __FUNCT__ 2121 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2122 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2123 { 2124 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2125 PetscErrorCode ierr; 2126 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2127 PetscScalar dx,*y,*xx,*w3_array; 2128 PetscScalar *vscale_array; 2129 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2130 Vec w1,w2,w3; 2131 void *fctx = coloring->fctx; 2132 PetscTruth flg = PETSC_FALSE; 2133 2134 PetscFunctionBegin; 2135 if (!coloring->w1) { 2136 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2137 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2138 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2139 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2140 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2141 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2142 } 2143 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2144 2145 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2146 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2147 if (flg) { 2148 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2149 } else { 2150 PetscTruth assembled; 2151 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2152 if (assembled) { 2153 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2154 } 2155 } 2156 2157 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2158 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2159 2160 /* 2161 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2162 coloring->F for the coarser grids from the finest 2163 */ 2164 if (coloring->F) { 2165 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2166 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2167 if (m1 != m2) { 2168 coloring->F = 0; 2169 } 2170 } 2171 2172 if (coloring->F) { 2173 w1 = coloring->F; 2174 coloring->F = 0; 2175 } else { 2176 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2177 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2178 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2179 } 2180 2181 /* 2182 Compute all the scale factors and share with other processors 2183 */ 2184 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2185 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2186 for (k=0; k<coloring->ncolors; k++) { 2187 /* 2188 Loop over each column associated with color adding the 2189 perturbation to the vector w3. 2190 */ 2191 for (l=0; l<coloring->ncolumns[k]; l++) { 2192 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2193 dx = xx[col]; 2194 if (dx == 0.0) dx = 1.0; 2195 #if !defined(PETSC_USE_COMPLEX) 2196 if (dx < umin && dx >= 0.0) dx = umin; 2197 else if (dx < 0.0 && dx > -umin) dx = -umin; 2198 #else 2199 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2200 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2201 #endif 2202 dx *= epsilon; 2203 vscale_array[col] = 1.0/dx; 2204 } 2205 } 2206 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2207 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2208 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2209 2210 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2211 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2212 2213 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2214 else vscaleforrow = coloring->columnsforrow; 2215 2216 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2217 /* 2218 Loop over each color 2219 */ 2220 for (k=0; k<coloring->ncolors; k++) { 2221 coloring->currentcolor = k; 2222 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2223 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2224 /* 2225 Loop over each column associated with color adding the 2226 perturbation to the vector w3. 2227 */ 2228 for (l=0; l<coloring->ncolumns[k]; l++) { 2229 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2230 dx = xx[col]; 2231 if (dx == 0.0) dx = 1.0; 2232 #if !defined(PETSC_USE_COMPLEX) 2233 if (dx < umin && dx >= 0.0) dx = umin; 2234 else if (dx < 0.0 && dx > -umin) dx = -umin; 2235 #else 2236 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2237 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2238 #endif 2239 dx *= epsilon; 2240 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2241 w3_array[col] += dx; 2242 } 2243 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2244 2245 /* 2246 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2247 */ 2248 2249 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2250 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2251 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2252 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2253 2254 /* 2255 Loop over rows of vector, putting results into Jacobian matrix 2256 */ 2257 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2258 for (l=0; l<coloring->nrows[k]; l++) { 2259 row = coloring->rows[k][l]; 2260 col = coloring->columnsforrow[k][l]; 2261 y[row] *= vscale_array[vscaleforrow[k][l]]; 2262 srow = row + start; 2263 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2264 } 2265 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2266 } 2267 coloring->currentcolor = k; 2268 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2269 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2270 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2271 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2272 PetscFunctionReturn(0); 2273 } 2274 2275 #undef __FUNCT__ 2276 #define __FUNCT__ "MatAXPY_SeqAIJ" 2277 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2278 { 2279 PetscErrorCode ierr; 2280 PetscInt i; 2281 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2282 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2283 2284 PetscFunctionBegin; 2285 if (str == SAME_NONZERO_PATTERN) { 2286 PetscScalar alpha = a; 2287 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2288 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2289 if (y->xtoy && y->XtoY != X) { 2290 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2291 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2292 } 2293 if (!y->xtoy) { /* get xtoy */ 2294 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2295 y->XtoY = X; 2296 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2297 } 2298 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2299 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2300 } else { 2301 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2302 } 2303 PetscFunctionReturn(0); 2304 } 2305 2306 #undef __FUNCT__ 2307 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2308 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2309 { 2310 PetscErrorCode ierr; 2311 2312 PetscFunctionBegin; 2313 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2314 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2315 PetscFunctionReturn(0); 2316 } 2317 2318 #undef __FUNCT__ 2319 #define __FUNCT__ "MatConjugate_SeqAIJ" 2320 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2321 { 2322 #if defined(PETSC_USE_COMPLEX) 2323 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2324 PetscInt i,nz; 2325 PetscScalar *a; 2326 2327 PetscFunctionBegin; 2328 nz = aij->nz; 2329 a = aij->a; 2330 for (i=0; i<nz; i++) { 2331 a[i] = PetscConj(a[i]); 2332 } 2333 #else 2334 PetscFunctionBegin; 2335 #endif 2336 PetscFunctionReturn(0); 2337 } 2338 2339 #undef __FUNCT__ 2340 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2341 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2342 { 2343 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2344 PetscErrorCode ierr; 2345 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2346 PetscReal atmp; 2347 PetscScalar *x; 2348 MatScalar *aa; 2349 2350 PetscFunctionBegin; 2351 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2352 aa = a->a; 2353 ai = a->i; 2354 aj = a->j; 2355 2356 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2357 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2358 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2359 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2360 for (i=0; i<m; i++) { 2361 ncols = ai[1] - ai[0]; ai++; 2362 x[i] = 0.0; 2363 for (j=0; j<ncols; j++){ 2364 atmp = PetscAbsScalar(*aa); 2365 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2366 aa++; aj++; 2367 } 2368 } 2369 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2370 PetscFunctionReturn(0); 2371 } 2372 2373 #undef __FUNCT__ 2374 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2375 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2376 { 2377 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2378 PetscErrorCode ierr; 2379 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2380 PetscScalar *x; 2381 MatScalar *aa; 2382 2383 PetscFunctionBegin; 2384 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2385 aa = a->a; 2386 ai = a->i; 2387 aj = a->j; 2388 2389 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2390 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2391 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2392 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2393 for (i=0; i<m; i++) { 2394 ncols = ai[1] - ai[0]; ai++; 2395 if (ncols == A->cmap->n) { /* row is dense */ 2396 x[i] = *aa; if (idx) idx[i] = 0; 2397 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2398 x[i] = 0.0; 2399 if (idx) { 2400 idx[i] = 0; /* in case ncols is zero */ 2401 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2402 if (aj[j] > j) { 2403 idx[i] = j; 2404 break; 2405 } 2406 } 2407 } 2408 } 2409 for (j=0; j<ncols; j++){ 2410 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2411 aa++; aj++; 2412 } 2413 } 2414 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2415 PetscFunctionReturn(0); 2416 } 2417 2418 #undef __FUNCT__ 2419 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2420 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2421 { 2422 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2423 PetscErrorCode ierr; 2424 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2425 PetscReal atmp; 2426 PetscScalar *x; 2427 MatScalar *aa; 2428 2429 PetscFunctionBegin; 2430 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2431 aa = a->a; 2432 ai = a->i; 2433 aj = a->j; 2434 2435 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2436 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2437 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2438 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2439 for (i=0; i<m; i++) { 2440 ncols = ai[1] - ai[0]; ai++; 2441 if (ncols) { 2442 /* Get first nonzero */ 2443 for(j = 0; j < ncols; j++) { 2444 atmp = PetscAbsScalar(aa[j]); 2445 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2446 } 2447 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2448 } else { 2449 x[i] = 0.0; if (idx) idx[i] = 0; 2450 } 2451 for(j = 0; j < ncols; j++) { 2452 atmp = PetscAbsScalar(*aa); 2453 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2454 aa++; aj++; 2455 } 2456 } 2457 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2458 PetscFunctionReturn(0); 2459 } 2460 2461 #undef __FUNCT__ 2462 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2463 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2464 { 2465 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2466 PetscErrorCode ierr; 2467 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2468 PetscScalar *x; 2469 MatScalar *aa; 2470 2471 PetscFunctionBegin; 2472 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2473 aa = a->a; 2474 ai = a->i; 2475 aj = a->j; 2476 2477 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2478 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2479 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2480 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2481 for (i=0; i<m; i++) { 2482 ncols = ai[1] - ai[0]; ai++; 2483 if (ncols == A->cmap->n) { /* row is dense */ 2484 x[i] = *aa; if (idx) idx[i] = 0; 2485 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2486 x[i] = 0.0; 2487 if (idx) { /* find first implicit 0.0 in the row */ 2488 idx[i] = 0; /* in case ncols is zero */ 2489 for (j=0;j<ncols;j++) { 2490 if (aj[j] > j) { 2491 idx[i] = j; 2492 break; 2493 } 2494 } 2495 } 2496 } 2497 for (j=0; j<ncols; j++){ 2498 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2499 aa++; aj++; 2500 } 2501 } 2502 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2503 PetscFunctionReturn(0); 2504 } 2505 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2506 /* -------------------------------------------------------------------*/ 2507 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2508 MatGetRow_SeqAIJ, 2509 MatRestoreRow_SeqAIJ, 2510 MatMult_SeqAIJ, 2511 /* 4*/ MatMultAdd_SeqAIJ, 2512 MatMultTranspose_SeqAIJ, 2513 MatMultTransposeAdd_SeqAIJ, 2514 0, 2515 0, 2516 0, 2517 /*10*/ 0, 2518 MatLUFactor_SeqAIJ, 2519 0, 2520 MatSOR_SeqAIJ, 2521 MatTranspose_SeqAIJ, 2522 /*15*/ MatGetInfo_SeqAIJ, 2523 MatEqual_SeqAIJ, 2524 MatGetDiagonal_SeqAIJ, 2525 MatDiagonalScale_SeqAIJ, 2526 MatNorm_SeqAIJ, 2527 /*20*/ 0, 2528 MatAssemblyEnd_SeqAIJ, 2529 MatSetOption_SeqAIJ, 2530 MatZeroEntries_SeqAIJ, 2531 /*24*/ MatZeroRows_SeqAIJ, 2532 0, 2533 0, 2534 0, 2535 0, 2536 /*29*/ MatSetUpPreallocation_SeqAIJ, 2537 0, 2538 0, 2539 MatGetArray_SeqAIJ, 2540 MatRestoreArray_SeqAIJ, 2541 /*34*/ MatDuplicate_SeqAIJ, 2542 0, 2543 0, 2544 MatILUFactor_SeqAIJ, 2545 0, 2546 /*39*/ MatAXPY_SeqAIJ, 2547 MatGetSubMatrices_SeqAIJ, 2548 MatIncreaseOverlap_SeqAIJ, 2549 MatGetValues_SeqAIJ, 2550 MatCopy_SeqAIJ, 2551 /*44*/ MatGetRowMax_SeqAIJ, 2552 MatScale_SeqAIJ, 2553 0, 2554 MatDiagonalSet_SeqAIJ, 2555 0, 2556 /*49*/ MatSetBlockSize_SeqAIJ, 2557 MatGetRowIJ_SeqAIJ, 2558 MatRestoreRowIJ_SeqAIJ, 2559 MatGetColumnIJ_SeqAIJ, 2560 MatRestoreColumnIJ_SeqAIJ, 2561 /*54*/ MatFDColoringCreate_SeqAIJ, 2562 0, 2563 0, 2564 MatPermute_SeqAIJ, 2565 0, 2566 /*59*/ 0, 2567 MatDestroy_SeqAIJ, 2568 MatView_SeqAIJ, 2569 0, 2570 0, 2571 /*64*/ 0, 2572 0, 2573 0, 2574 0, 2575 0, 2576 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2577 MatGetRowMinAbs_SeqAIJ, 2578 0, 2579 MatSetColoring_SeqAIJ, 2580 #if defined(PETSC_HAVE_ADIC) 2581 MatSetValuesAdic_SeqAIJ, 2582 #else 2583 0, 2584 #endif 2585 /*74*/ MatSetValuesAdifor_SeqAIJ, 2586 MatFDColoringApply_AIJ, 2587 0, 2588 0, 2589 0, 2590 /*79*/ 0, 2591 0, 2592 0, 2593 0, 2594 MatLoad_SeqAIJ, 2595 /*84*/ MatIsSymmetric_SeqAIJ, 2596 MatIsHermitian_SeqAIJ, 2597 0, 2598 0, 2599 0, 2600 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2601 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2602 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2603 MatPtAP_Basic, 2604 MatPtAPSymbolic_SeqAIJ, 2605 /*94*/ MatPtAPNumeric_SeqAIJ, 2606 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2607 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2608 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2609 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2610 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2611 0, 2612 0, 2613 MatConjugate_SeqAIJ, 2614 0, 2615 /*104*/MatSetValuesRow_SeqAIJ, 2616 MatRealPart_SeqAIJ, 2617 MatImaginaryPart_SeqAIJ, 2618 0, 2619 0, 2620 /*109*/0, 2621 0, 2622 MatGetRowMin_SeqAIJ, 2623 0, 2624 MatMissingDiagonal_SeqAIJ, 2625 /*114*/0, 2626 0, 2627 0, 2628 0, 2629 0, 2630 /*119*/0, 2631 0, 2632 0, 2633 0 2634 }; 2635 2636 EXTERN_C_BEGIN 2637 #undef __FUNCT__ 2638 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2639 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2640 { 2641 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2642 PetscInt i,nz,n; 2643 2644 PetscFunctionBegin; 2645 2646 nz = aij->maxnz; 2647 n = mat->rmap->n; 2648 for (i=0; i<nz; i++) { 2649 aij->j[i] = indices[i]; 2650 } 2651 aij->nz = nz; 2652 for (i=0; i<n; i++) { 2653 aij->ilen[i] = aij->imax[i]; 2654 } 2655 2656 PetscFunctionReturn(0); 2657 } 2658 EXTERN_C_END 2659 2660 #undef __FUNCT__ 2661 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2662 /*@ 2663 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2664 in the matrix. 2665 2666 Input Parameters: 2667 + mat - the SeqAIJ matrix 2668 - indices - the column indices 2669 2670 Level: advanced 2671 2672 Notes: 2673 This can be called if you have precomputed the nonzero structure of the 2674 matrix and want to provide it to the matrix object to improve the performance 2675 of the MatSetValues() operation. 2676 2677 You MUST have set the correct numbers of nonzeros per row in the call to 2678 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2679 2680 MUST be called before any calls to MatSetValues(); 2681 2682 The indices should start with zero, not one. 2683 2684 @*/ 2685 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2686 { 2687 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2688 2689 PetscFunctionBegin; 2690 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2691 PetscValidPointer(indices,2); 2692 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2693 if (f) { 2694 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2695 } else { 2696 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2697 } 2698 PetscFunctionReturn(0); 2699 } 2700 2701 /* ----------------------------------------------------------------------------------------*/ 2702 2703 EXTERN_C_BEGIN 2704 #undef __FUNCT__ 2705 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2706 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2707 { 2708 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2709 PetscErrorCode ierr; 2710 size_t nz = aij->i[mat->rmap->n]; 2711 2712 PetscFunctionBegin; 2713 if (aij->nonew != 1) { 2714 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2715 } 2716 2717 /* allocate space for values if not already there */ 2718 if (!aij->saved_values) { 2719 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2720 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2721 } 2722 2723 /* copy values over */ 2724 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2725 PetscFunctionReturn(0); 2726 } 2727 EXTERN_C_END 2728 2729 #undef __FUNCT__ 2730 #define __FUNCT__ "MatStoreValues" 2731 /*@ 2732 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2733 example, reuse of the linear part of a Jacobian, while recomputing the 2734 nonlinear portion. 2735 2736 Collect on Mat 2737 2738 Input Parameters: 2739 . mat - the matrix (currently only AIJ matrices support this option) 2740 2741 Level: advanced 2742 2743 Common Usage, with SNESSolve(): 2744 $ Create Jacobian matrix 2745 $ Set linear terms into matrix 2746 $ Apply boundary conditions to matrix, at this time matrix must have 2747 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2748 $ boundary conditions again will not change the nonzero structure 2749 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2750 $ ierr = MatStoreValues(mat); 2751 $ Call SNESSetJacobian() with matrix 2752 $ In your Jacobian routine 2753 $ ierr = MatRetrieveValues(mat); 2754 $ Set nonlinear terms in matrix 2755 2756 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2757 $ // build linear portion of Jacobian 2758 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2759 $ ierr = MatStoreValues(mat); 2760 $ loop over nonlinear iterations 2761 $ ierr = MatRetrieveValues(mat); 2762 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2763 $ // call MatAssemblyBegin/End() on matrix 2764 $ Solve linear system with Jacobian 2765 $ endloop 2766 2767 Notes: 2768 Matrix must already be assemblied before calling this routine 2769 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2770 calling this routine. 2771 2772 When this is called multiple times it overwrites the previous set of stored values 2773 and does not allocated additional space. 2774 2775 .seealso: MatRetrieveValues() 2776 2777 @*/ 2778 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2779 { 2780 PetscErrorCode ierr,(*f)(Mat); 2781 2782 PetscFunctionBegin; 2783 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2784 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2785 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2786 2787 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2788 if (f) { 2789 ierr = (*f)(mat);CHKERRQ(ierr); 2790 } else { 2791 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2792 } 2793 PetscFunctionReturn(0); 2794 } 2795 2796 EXTERN_C_BEGIN 2797 #undef __FUNCT__ 2798 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2799 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2800 { 2801 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2802 PetscErrorCode ierr; 2803 PetscInt nz = aij->i[mat->rmap->n]; 2804 2805 PetscFunctionBegin; 2806 if (aij->nonew != 1) { 2807 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2808 } 2809 if (!aij->saved_values) { 2810 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2811 } 2812 /* copy values over */ 2813 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2814 PetscFunctionReturn(0); 2815 } 2816 EXTERN_C_END 2817 2818 #undef __FUNCT__ 2819 #define __FUNCT__ "MatRetrieveValues" 2820 /*@ 2821 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2822 example, reuse of the linear part of a Jacobian, while recomputing the 2823 nonlinear portion. 2824 2825 Collect on Mat 2826 2827 Input Parameters: 2828 . mat - the matrix (currently on AIJ matrices support this option) 2829 2830 Level: advanced 2831 2832 .seealso: MatStoreValues() 2833 2834 @*/ 2835 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2836 { 2837 PetscErrorCode ierr,(*f)(Mat); 2838 2839 PetscFunctionBegin; 2840 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2841 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2842 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2843 2844 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2845 if (f) { 2846 ierr = (*f)(mat);CHKERRQ(ierr); 2847 } else { 2848 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2849 } 2850 PetscFunctionReturn(0); 2851 } 2852 2853 2854 /* --------------------------------------------------------------------------------*/ 2855 #undef __FUNCT__ 2856 #define __FUNCT__ "MatCreateSeqAIJ" 2857 /*@C 2858 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2859 (the default parallel PETSc format). For good matrix assembly performance 2860 the user should preallocate the matrix storage by setting the parameter nz 2861 (or the array nnz). By setting these parameters accurately, performance 2862 during matrix assembly can be increased by more than a factor of 50. 2863 2864 Collective on MPI_Comm 2865 2866 Input Parameters: 2867 + comm - MPI communicator, set to PETSC_COMM_SELF 2868 . m - number of rows 2869 . n - number of columns 2870 . nz - number of nonzeros per row (same for all rows) 2871 - nnz - array containing the number of nonzeros in the various rows 2872 (possibly different for each row) or PETSC_NULL 2873 2874 Output Parameter: 2875 . A - the matrix 2876 2877 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2878 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2879 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2880 2881 Notes: 2882 If nnz is given then nz is ignored 2883 2884 The AIJ format (also called the Yale sparse matrix format or 2885 compressed row storage), is fully compatible with standard Fortran 77 2886 storage. That is, the stored row and column indices can begin at 2887 either one (as in Fortran) or zero. See the users' manual for details. 2888 2889 Specify the preallocated storage with either nz or nnz (not both). 2890 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2891 allocation. For large problems you MUST preallocate memory or you 2892 will get TERRIBLE performance, see the users' manual chapter on matrices. 2893 2894 By default, this format uses inodes (identical nodes) when possible, to 2895 improve numerical efficiency of matrix-vector products and solves. We 2896 search for consecutive rows with the same nonzero structure, thereby 2897 reusing matrix information to achieve increased efficiency. 2898 2899 Options Database Keys: 2900 + -mat_no_inode - Do not use inodes 2901 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2902 2903 Level: intermediate 2904 2905 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2906 2907 @*/ 2908 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2909 { 2910 PetscErrorCode ierr; 2911 2912 PetscFunctionBegin; 2913 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2914 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2915 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2916 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2917 PetscFunctionReturn(0); 2918 } 2919 2920 #undef __FUNCT__ 2921 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2922 /*@C 2923 MatSeqAIJSetPreallocation - For good matrix assembly performance 2924 the user should preallocate the matrix storage by setting the parameter nz 2925 (or the array nnz). By setting these parameters accurately, performance 2926 during matrix assembly can be increased by more than a factor of 50. 2927 2928 Collective on MPI_Comm 2929 2930 Input Parameters: 2931 + B - The matrix-free 2932 . nz - number of nonzeros per row (same for all rows) 2933 - nnz - array containing the number of nonzeros in the various rows 2934 (possibly different for each row) or PETSC_NULL 2935 2936 Notes: 2937 If nnz is given then nz is ignored 2938 2939 The AIJ format (also called the Yale sparse matrix format or 2940 compressed row storage), is fully compatible with standard Fortran 77 2941 storage. That is, the stored row and column indices can begin at 2942 either one (as in Fortran) or zero. See the users' manual for details. 2943 2944 Specify the preallocated storage with either nz or nnz (not both). 2945 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2946 allocation. For large problems you MUST preallocate memory or you 2947 will get TERRIBLE performance, see the users' manual chapter on matrices. 2948 2949 You can call MatGetInfo() to get information on how effective the preallocation was; 2950 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2951 You can also run with the option -info and look for messages with the string 2952 malloc in them to see if additional memory allocation was needed. 2953 2954 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2955 entries or columns indices 2956 2957 By default, this format uses inodes (identical nodes) when possible, to 2958 improve numerical efficiency of matrix-vector products and solves. We 2959 search for consecutive rows with the same nonzero structure, thereby 2960 reusing matrix information to achieve increased efficiency. 2961 2962 Options Database Keys: 2963 + -mat_no_inode - Do not use inodes 2964 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2965 - -mat_aij_oneindex - Internally use indexing starting at 1 2966 rather than 0. Note that when calling MatSetValues(), 2967 the user still MUST index entries starting at 0! 2968 2969 Level: intermediate 2970 2971 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2972 2973 @*/ 2974 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2975 { 2976 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2977 2978 PetscFunctionBegin; 2979 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2980 if (f) { 2981 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2982 } 2983 PetscFunctionReturn(0); 2984 } 2985 2986 EXTERN_C_BEGIN 2987 #undef __FUNCT__ 2988 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2989 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2990 { 2991 Mat_SeqAIJ *b; 2992 PetscTruth skipallocation = PETSC_FALSE; 2993 PetscErrorCode ierr; 2994 PetscInt i; 2995 2996 PetscFunctionBegin; 2997 2998 if (nz == MAT_SKIP_ALLOCATION) { 2999 skipallocation = PETSC_TRUE; 3000 nz = 0; 3001 } 3002 3003 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3004 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3005 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3006 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3007 3008 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3009 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3010 if (nnz) { 3011 for (i=0; i<B->rmap->n; i++) { 3012 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]); 3013 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); 3014 } 3015 } 3016 3017 B->preallocated = PETSC_TRUE; 3018 b = (Mat_SeqAIJ*)B->data; 3019 3020 if (!skipallocation) { 3021 if (!b->imax) { 3022 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3023 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3024 } 3025 if (!nnz) { 3026 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3027 else if (nz <= 0) nz = 1; 3028 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3029 nz = nz*B->rmap->n; 3030 } else { 3031 nz = 0; 3032 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3033 } 3034 /* b->ilen will count nonzeros in each row so far. */ 3035 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3036 3037 /* allocate the matrix space */ 3038 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3039 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3040 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3041 b->i[0] = 0; 3042 for (i=1; i<B->rmap->n+1; i++) { 3043 b->i[i] = b->i[i-1] + b->imax[i-1]; 3044 } 3045 b->singlemalloc = PETSC_TRUE; 3046 b->free_a = PETSC_TRUE; 3047 b->free_ij = PETSC_TRUE; 3048 } else { 3049 b->free_a = PETSC_FALSE; 3050 b->free_ij = PETSC_FALSE; 3051 } 3052 3053 b->nz = 0; 3054 b->maxnz = nz; 3055 B->info.nz_unneeded = (double)b->maxnz; 3056 PetscFunctionReturn(0); 3057 } 3058 EXTERN_C_END 3059 3060 #undef __FUNCT__ 3061 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3062 /*@ 3063 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3064 3065 Input Parameters: 3066 + B - the matrix 3067 . i - the indices into j for the start of each row (starts with zero) 3068 . j - the column indices for each row (starts with zero) these must be sorted for each row 3069 - v - optional values in the matrix 3070 3071 Level: developer 3072 3073 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3074 3075 .keywords: matrix, aij, compressed row, sparse, sequential 3076 3077 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3078 @*/ 3079 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3080 { 3081 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3082 PetscErrorCode ierr; 3083 3084 PetscFunctionBegin; 3085 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3086 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3087 if (f) { 3088 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3089 } 3090 PetscFunctionReturn(0); 3091 } 3092 3093 EXTERN_C_BEGIN 3094 #undef __FUNCT__ 3095 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3096 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3097 { 3098 PetscInt i; 3099 PetscInt m,n; 3100 PetscInt nz; 3101 PetscInt *nnz, nz_max = 0; 3102 PetscScalar *values; 3103 PetscErrorCode ierr; 3104 3105 PetscFunctionBegin; 3106 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3107 3108 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3109 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3110 for(i = 0; i < m; i++) { 3111 nz = Ii[i+1]- Ii[i]; 3112 nz_max = PetscMax(nz_max, nz); 3113 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3114 nnz[i] = nz; 3115 } 3116 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3117 ierr = PetscFree(nnz);CHKERRQ(ierr); 3118 3119 if (v) { 3120 values = (PetscScalar*) v; 3121 } else { 3122 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3123 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3124 } 3125 3126 for(i = 0; i < m; i++) { 3127 nz = Ii[i+1] - Ii[i]; 3128 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3129 } 3130 3131 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3132 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3133 3134 if (!v) { 3135 ierr = PetscFree(values);CHKERRQ(ierr); 3136 } 3137 PetscFunctionReturn(0); 3138 } 3139 EXTERN_C_END 3140 3141 #include "../src/mat/impls/dense/seq/dense.h" 3142 #include "private/petscaxpy.h" 3143 3144 #undef __FUNCT__ 3145 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3146 /* 3147 Computes (B'*A')' since computing B*A directly is untenable 3148 3149 n p p 3150 ( ) ( ) ( ) 3151 m ( A ) * n ( B ) = m ( C ) 3152 ( ) ( ) ( ) 3153 3154 */ 3155 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3156 { 3157 PetscErrorCode ierr; 3158 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3159 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3160 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3161 PetscInt i,n,m,q,p; 3162 const PetscInt *ii,*idx; 3163 const PetscScalar *b,*a,*a_q; 3164 PetscScalar *c,*c_q; 3165 3166 PetscFunctionBegin; 3167 m = A->rmap->n; 3168 n = A->cmap->n; 3169 p = B->cmap->n; 3170 a = sub_a->v; 3171 b = sub_b->a; 3172 c = sub_c->v; 3173 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3174 3175 ii = sub_b->i; 3176 idx = sub_b->j; 3177 for (i=0; i<n; i++) { 3178 q = ii[i+1] - ii[i]; 3179 while (q-->0) { 3180 c_q = c + m*(*idx); 3181 a_q = a + m*i; 3182 PetscAXPY(c_q,*b,a_q,m); 3183 idx++; 3184 b++; 3185 } 3186 } 3187 PetscFunctionReturn(0); 3188 } 3189 3190 #undef __FUNCT__ 3191 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3192 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3193 { 3194 PetscErrorCode ierr; 3195 PetscInt m=A->rmap->n,n=B->cmap->n; 3196 Mat Cmat; 3197 3198 PetscFunctionBegin; 3199 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); 3200 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3201 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3202 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3203 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3204 Cmat->assembled = PETSC_TRUE; 3205 *C = Cmat; 3206 PetscFunctionReturn(0); 3207 } 3208 3209 /* ----------------------------------------------------------------*/ 3210 #undef __FUNCT__ 3211 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3212 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3213 { 3214 PetscErrorCode ierr; 3215 3216 PetscFunctionBegin; 3217 if (scall == MAT_INITIAL_MATRIX){ 3218 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3219 } 3220 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3221 PetscFunctionReturn(0); 3222 } 3223 3224 3225 /*MC 3226 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3227 based on compressed sparse row format. 3228 3229 Options Database Keys: 3230 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3231 3232 Level: beginner 3233 3234 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3235 M*/ 3236 3237 EXTERN_C_BEGIN 3238 #if defined(PETSC_HAVE_PASTIX) 3239 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3240 #endif 3241 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3242 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3243 #endif 3244 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3246 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3247 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3248 #if defined(PETSC_HAVE_MUMPS) 3249 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3250 #endif 3251 #if defined(PETSC_HAVE_SUPERLU) 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3253 #endif 3254 #if defined(PETSC_HAVE_SUPERLU_DIST) 3255 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3256 #endif 3257 #if defined(PETSC_HAVE_SPOOLES) 3258 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3259 #endif 3260 #if defined(PETSC_HAVE_UMFPACK) 3261 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3262 #endif 3263 #if defined(PETSC_HAVE_CHOLMOD) 3264 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3265 #endif 3266 #if defined(PETSC_HAVE_LUSOL) 3267 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3268 #endif 3269 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3270 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3271 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3272 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3273 #endif 3274 EXTERN_C_END 3275 3276 3277 EXTERN_C_BEGIN 3278 #undef __FUNCT__ 3279 #define __FUNCT__ "MatCreate_SeqAIJ" 3280 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3281 { 3282 Mat_SeqAIJ *b; 3283 PetscErrorCode ierr; 3284 PetscMPIInt size; 3285 3286 PetscFunctionBegin; 3287 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3288 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3289 3290 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3291 B->data = (void*)b; 3292 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3293 B->mapping = 0; 3294 b->row = 0; 3295 b->col = 0; 3296 b->icol = 0; 3297 b->reallocs = 0; 3298 b->ignorezeroentries = PETSC_FALSE; 3299 b->roworiented = PETSC_TRUE; 3300 b->nonew = 0; 3301 b->diag = 0; 3302 b->solve_work = 0; 3303 B->spptr = 0; 3304 b->saved_values = 0; 3305 b->idiag = 0; 3306 b->mdiag = 0; 3307 b->ssor_work = 0; 3308 b->omega = 1.0; 3309 b->fshift = 0.0; 3310 b->idiagvalid = PETSC_FALSE; 3311 b->keepnonzeropattern = PETSC_FALSE; 3312 b->xtoy = 0; 3313 b->XtoY = 0; 3314 b->compressedrow.use = PETSC_FALSE; 3315 b->compressedrow.nrows = B->rmap->n; 3316 b->compressedrow.i = PETSC_NULL; 3317 b->compressedrow.rindex = PETSC_NULL; 3318 b->compressedrow.checked = PETSC_FALSE; 3319 B->same_nonzero = PETSC_FALSE; 3320 3321 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3322 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3324 "MatGetFactor_seqaij_matlab", 3325 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3326 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3328 #endif 3329 #if defined(PETSC_HAVE_PASTIX) 3330 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3331 "MatGetFactor_seqaij_pastix", 3332 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3333 #endif 3334 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3336 "MatGetFactor_seqaij_essl", 3337 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3338 #endif 3339 #if defined(PETSC_HAVE_SUPERLU) 3340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3341 "MatGetFactor_seqaij_superlu", 3342 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3343 #endif 3344 #if defined(PETSC_HAVE_SUPERLU_DIST) 3345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3346 "MatGetFactor_seqaij_superlu_dist", 3347 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3348 #endif 3349 #if defined(PETSC_HAVE_SPOOLES) 3350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3351 "MatGetFactor_seqaij_spooles", 3352 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3353 #endif 3354 #if defined(PETSC_HAVE_MUMPS) 3355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3356 "MatGetFactor_aij_mumps", 3357 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3358 #endif 3359 #if defined(PETSC_HAVE_UMFPACK) 3360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3361 "MatGetFactor_seqaij_umfpack", 3362 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3363 #endif 3364 #if defined(PETSC_HAVE_CHOLMOD) 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3366 "MatGetFactor_seqaij_cholmod", 3367 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3368 #endif 3369 #if defined(PETSC_HAVE_LUSOL) 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3371 "MatGetFactor_seqaij_lusol", 3372 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3373 #endif 3374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3375 "MatGetFactor_seqaij_petsc", 3376 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3377 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3378 "MatGetFactorAvailable_seqaij_petsc", 3379 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3380 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3381 "MatGetFactor_seqaij_bas", 3382 MatGetFactor_seqaij_bas); 3383 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3384 "MatSeqAIJSetColumnIndices_SeqAIJ", 3385 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3387 "MatStoreValues_SeqAIJ", 3388 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3390 "MatRetrieveValues_SeqAIJ", 3391 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3392 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3393 "MatConvert_SeqAIJ_SeqSBAIJ", 3394 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3395 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3396 "MatConvert_SeqAIJ_SeqBAIJ", 3397 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3398 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C", 3399 "MatConvert_SeqAIJ_SeqAIJPERM", 3400 MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3401 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C", 3402 "MatConvert_SeqAIJ_SeqAIJCRL", 3403 MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3404 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3405 "MatIsTranspose_SeqAIJ", 3406 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3407 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3408 "MatIsHermitianTranspose_SeqAIJ", 3409 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3410 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3411 "MatSeqAIJSetPreallocation_SeqAIJ", 3412 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3413 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3414 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3415 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3416 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3417 "MatReorderForNonzeroDiagonal_SeqAIJ", 3418 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3419 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3420 "MatMatMult_SeqDense_SeqAIJ", 3421 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3422 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3423 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3424 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3425 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3426 "MatMatMultNumeric_SeqDense_SeqAIJ", 3427 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3428 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3429 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3430 PetscFunctionReturn(0); 3431 } 3432 EXTERN_C_END 3433 3434 #undef __FUNCT__ 3435 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3436 /* 3437 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3438 */ 3439 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3440 { 3441 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3442 PetscErrorCode ierr; 3443 PetscInt i,m = A->rmap->n; 3444 3445 PetscFunctionBegin; 3446 c = (Mat_SeqAIJ*)C->data; 3447 3448 C->factortype = A->factortype; 3449 c->row = 0; 3450 c->col = 0; 3451 c->icol = 0; 3452 c->reallocs = 0; 3453 3454 C->assembled = PETSC_TRUE; 3455 3456 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3457 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3458 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3459 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3460 3461 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3462 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3463 for (i=0; i<m; i++) { 3464 c->imax[i] = a->imax[i]; 3465 c->ilen[i] = a->ilen[i]; 3466 } 3467 3468 /* allocate the matrix space */ 3469 if (mallocmatspace){ 3470 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3471 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3472 c->singlemalloc = PETSC_TRUE; 3473 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3474 if (m > 0) { 3475 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3476 if (cpvalues == MAT_COPY_VALUES) { 3477 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3478 } else { 3479 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3480 } 3481 } 3482 } 3483 3484 c->ignorezeroentries = a->ignorezeroentries; 3485 c->roworiented = a->roworiented; 3486 c->nonew = a->nonew; 3487 if (a->diag) { 3488 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3489 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3490 for (i=0; i<m; i++) { 3491 c->diag[i] = a->diag[i]; 3492 } 3493 } else c->diag = 0; 3494 c->solve_work = 0; 3495 c->saved_values = 0; 3496 c->idiag = 0; 3497 c->ssor_work = 0; 3498 c->keepnonzeropattern = a->keepnonzeropattern; 3499 c->free_a = PETSC_TRUE; 3500 c->free_ij = PETSC_TRUE; 3501 c->xtoy = 0; 3502 c->XtoY = 0; 3503 3504 c->nz = a->nz; 3505 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3506 C->preallocated = PETSC_TRUE; 3507 3508 c->compressedrow.use = a->compressedrow.use; 3509 c->compressedrow.nrows = a->compressedrow.nrows; 3510 c->compressedrow.checked = a->compressedrow.checked; 3511 if (a->compressedrow.checked && a->compressedrow.use){ 3512 i = a->compressedrow.nrows; 3513 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3514 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3515 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3516 } else { 3517 c->compressedrow.use = PETSC_FALSE; 3518 c->compressedrow.i = PETSC_NULL; 3519 c->compressedrow.rindex = PETSC_NULL; 3520 } 3521 C->same_nonzero = A->same_nonzero; 3522 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3523 3524 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3525 PetscFunctionReturn(0); 3526 } 3527 3528 #undef __FUNCT__ 3529 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3530 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3531 { 3532 PetscErrorCode ierr; 3533 3534 PetscFunctionBegin; 3535 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3536 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3537 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3538 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3539 PetscFunctionReturn(0); 3540 } 3541 3542 #undef __FUNCT__ 3543 #define __FUNCT__ "MatLoad_SeqAIJ" 3544 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3545 { 3546 Mat_SeqAIJ *a; 3547 PetscErrorCode ierr; 3548 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3549 int fd; 3550 PetscMPIInt size; 3551 MPI_Comm comm; 3552 3553 PetscFunctionBegin; 3554 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3555 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3556 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3557 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3558 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3559 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3560 M = header[1]; N = header[2]; nz = header[3]; 3561 3562 if (nz < 0) { 3563 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3564 } 3565 3566 /* read in row lengths */ 3567 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3568 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3569 3570 /* check if sum of rowlengths is same as nz */ 3571 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3572 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); 3573 3574 /* set global size if not set already*/ 3575 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3576 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3577 } else { 3578 /* if sizes and type are already set, check if the vector global sizes are correct */ 3579 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3580 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); 3581 } 3582 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3583 a = (Mat_SeqAIJ*)newMat->data; 3584 3585 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3586 3587 /* read in nonzero values */ 3588 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3589 3590 /* set matrix "i" values */ 3591 a->i[0] = 0; 3592 for (i=1; i<= M; i++) { 3593 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3594 a->ilen[i-1] = rowlengths[i-1]; 3595 } 3596 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3597 3598 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3599 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3600 PetscFunctionReturn(0); 3601 } 3602 3603 #undef __FUNCT__ 3604 #define __FUNCT__ "MatEqual_SeqAIJ" 3605 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3606 { 3607 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3608 PetscErrorCode ierr; 3609 #if defined(PETSC_USE_COMPLEX) 3610 PetscInt k; 3611 #endif 3612 3613 PetscFunctionBegin; 3614 /* If the matrix dimensions are not equal,or no of nonzeros */ 3615 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3616 *flg = PETSC_FALSE; 3617 PetscFunctionReturn(0); 3618 } 3619 3620 /* if the a->i are the same */ 3621 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3622 if (!*flg) PetscFunctionReturn(0); 3623 3624 /* if a->j are the same */ 3625 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3626 if (!*flg) PetscFunctionReturn(0); 3627 3628 /* if a->a are the same */ 3629 #if defined(PETSC_USE_COMPLEX) 3630 for (k=0; k<a->nz; k++){ 3631 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3632 *flg = PETSC_FALSE; 3633 PetscFunctionReturn(0); 3634 } 3635 } 3636 #else 3637 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3638 #endif 3639 PetscFunctionReturn(0); 3640 } 3641 3642 #undef __FUNCT__ 3643 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3644 /*@ 3645 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3646 provided by the user. 3647 3648 Collective on MPI_Comm 3649 3650 Input Parameters: 3651 + comm - must be an MPI communicator of size 1 3652 . m - number of rows 3653 . n - number of columns 3654 . i - row indices 3655 . j - column indices 3656 - a - matrix values 3657 3658 Output Parameter: 3659 . mat - the matrix 3660 3661 Level: intermediate 3662 3663 Notes: 3664 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3665 once the matrix is destroyed 3666 3667 You cannot set new nonzero locations into this matrix, that will generate an error. 3668 3669 The i and j indices are 0 based 3670 3671 The format which is used for the sparse matrix input, is equivalent to a 3672 row-major ordering.. i.e for the following matrix, the input data expected is 3673 as shown: 3674 3675 1 0 0 3676 2 0 3 3677 4 5 6 3678 3679 i = {0,1,3,6} [size = nrow+1 = 3+1] 3680 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3681 v = {1,2,3,4,5,6} [size = nz = 6] 3682 3683 3684 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3685 3686 @*/ 3687 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3688 { 3689 PetscErrorCode ierr; 3690 PetscInt ii; 3691 Mat_SeqAIJ *aij; 3692 #if defined(PETSC_USE_DEBUG) 3693 PetscInt jj; 3694 #endif 3695 3696 PetscFunctionBegin; 3697 if (i[0]) { 3698 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3699 } 3700 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3701 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3702 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3703 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3704 aij = (Mat_SeqAIJ*)(*mat)->data; 3705 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3706 3707 aij->i = i; 3708 aij->j = j; 3709 aij->a = a; 3710 aij->singlemalloc = PETSC_FALSE; 3711 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3712 aij->free_a = PETSC_FALSE; 3713 aij->free_ij = PETSC_FALSE; 3714 3715 for (ii=0; ii<m; ii++) { 3716 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3717 #if defined(PETSC_USE_DEBUG) 3718 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]); 3719 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3720 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); 3721 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); 3722 } 3723 #endif 3724 } 3725 #if defined(PETSC_USE_DEBUG) 3726 for (ii=0; ii<aij->i[m]; ii++) { 3727 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3728 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]); 3729 } 3730 #endif 3731 3732 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3733 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3734 PetscFunctionReturn(0); 3735 } 3736 3737 #undef __FUNCT__ 3738 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3739 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3740 { 3741 PetscErrorCode ierr; 3742 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3743 3744 PetscFunctionBegin; 3745 if (coloring->ctype == IS_COLORING_GLOBAL) { 3746 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3747 a->coloring = coloring; 3748 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3749 PetscInt i,*larray; 3750 ISColoring ocoloring; 3751 ISColoringValue *colors; 3752 3753 /* set coloring for diagonal portion */ 3754 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3755 for (i=0; i<A->cmap->n; i++) { 3756 larray[i] = i; 3757 } 3758 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3759 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3760 for (i=0; i<A->cmap->n; i++) { 3761 colors[i] = coloring->colors[larray[i]]; 3762 } 3763 ierr = PetscFree(larray);CHKERRQ(ierr); 3764 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3765 a->coloring = ocoloring; 3766 } 3767 PetscFunctionReturn(0); 3768 } 3769 3770 #if defined(PETSC_HAVE_ADIC) 3771 EXTERN_C_BEGIN 3772 #include "adic/ad_utils.h" 3773 EXTERN_C_END 3774 3775 #undef __FUNCT__ 3776 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3777 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3778 { 3779 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3780 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3781 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3782 ISColoringValue *color; 3783 3784 PetscFunctionBegin; 3785 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3786 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3787 color = a->coloring->colors; 3788 /* loop over rows */ 3789 for (i=0; i<m; i++) { 3790 nz = ii[i+1] - ii[i]; 3791 /* loop over columns putting computed value into matrix */ 3792 for (j=0; j<nz; j++) { 3793 *v++ = values[color[*jj++]]; 3794 } 3795 values += nlen; /* jump to next row of derivatives */ 3796 } 3797 PetscFunctionReturn(0); 3798 } 3799 #endif 3800 3801 #undef __FUNCT__ 3802 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3803 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3804 { 3805 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3806 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3807 MatScalar *v = a->a; 3808 PetscScalar *values = (PetscScalar *)advalues; 3809 ISColoringValue *color; 3810 3811 PetscFunctionBegin; 3812 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3813 color = a->coloring->colors; 3814 /* loop over rows */ 3815 for (i=0; i<m; i++) { 3816 nz = ii[i+1] - ii[i]; 3817 /* loop over columns putting computed value into matrix */ 3818 for (j=0; j<nz; j++) { 3819 *v++ = values[color[*jj++]]; 3820 } 3821 values += nl; /* jump to next row of derivatives */ 3822 } 3823 PetscFunctionReturn(0); 3824 } 3825 3826 /* 3827 Special version for direct calls from Fortran 3828 */ 3829 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3830 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3831 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3832 #define matsetvaluesseqaij_ matsetvaluesseqaij 3833 #endif 3834 3835 /* Change these macros so can be used in void function */ 3836 #undef CHKERRQ 3837 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3838 #undef SETERRQ2 3839 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3840 3841 EXTERN_C_BEGIN 3842 #undef __FUNCT__ 3843 #define __FUNCT__ "matsetvaluesseqaij_" 3844 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3845 { 3846 Mat A = *AA; 3847 PetscInt m = *mm, n = *nn; 3848 InsertMode is = *isis; 3849 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3850 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3851 PetscInt *imax,*ai,*ailen; 3852 PetscErrorCode ierr; 3853 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3854 MatScalar *ap,value,*aa; 3855 PetscTruth ignorezeroentries = a->ignorezeroentries; 3856 PetscTruth roworiented = a->roworiented; 3857 3858 PetscFunctionBegin; 3859 ierr = MatPreallocated(A);CHKERRQ(ierr); 3860 imax = a->imax; 3861 ai = a->i; 3862 ailen = a->ilen; 3863 aj = a->j; 3864 aa = a->a; 3865 3866 for (k=0; k<m; k++) { /* loop over added rows */ 3867 row = im[k]; 3868 if (row < 0) continue; 3869 #if defined(PETSC_USE_DEBUG) 3870 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3871 #endif 3872 rp = aj + ai[row]; ap = aa + ai[row]; 3873 rmax = imax[row]; nrow = ailen[row]; 3874 low = 0; 3875 high = nrow; 3876 for (l=0; l<n; l++) { /* loop over added columns */ 3877 if (in[l] < 0) continue; 3878 #if defined(PETSC_USE_DEBUG) 3879 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3880 #endif 3881 col = in[l]; 3882 if (roworiented) { 3883 value = v[l + k*n]; 3884 } else { 3885 value = v[k + l*m]; 3886 } 3887 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3888 3889 if (col <= lastcol) low = 0; else high = nrow; 3890 lastcol = col; 3891 while (high-low > 5) { 3892 t = (low+high)/2; 3893 if (rp[t] > col) high = t; 3894 else low = t; 3895 } 3896 for (i=low; i<high; i++) { 3897 if (rp[i] > col) break; 3898 if (rp[i] == col) { 3899 if (is == ADD_VALUES) ap[i] += value; 3900 else ap[i] = value; 3901 goto noinsert; 3902 } 3903 } 3904 if (value == 0.0 && ignorezeroentries) goto noinsert; 3905 if (nonew == 1) goto noinsert; 3906 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3907 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3908 N = nrow++ - 1; a->nz++; high++; 3909 /* shift up all the later entries in this row */ 3910 for (ii=N; ii>=i; ii--) { 3911 rp[ii+1] = rp[ii]; 3912 ap[ii+1] = ap[ii]; 3913 } 3914 rp[i] = col; 3915 ap[i] = value; 3916 noinsert:; 3917 low = i + 1; 3918 } 3919 ailen[row] = nrow; 3920 } 3921 A->same_nonzero = PETSC_FALSE; 3922 PetscFunctionReturnVoid(); 3923 } 3924 EXTERN_C_END 3925