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