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