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]) { 2167 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2168 } 2169 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2170 } else { 2171 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2172 } 2173 PetscFunctionReturn(0); 2174 } 2175 2176 #undef __FUNCT__ 2177 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2178 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2179 { 2180 PetscErrorCode ierr; 2181 2182 PetscFunctionBegin; 2183 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2184 PetscFunctionReturn(0); 2185 } 2186 2187 #undef __FUNCT__ 2188 #define __FUNCT__ "MatGetArray_SeqAIJ" 2189 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2190 { 2191 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2192 PetscFunctionBegin; 2193 *array = a->a; 2194 PetscFunctionReturn(0); 2195 } 2196 2197 #undef __FUNCT__ 2198 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2199 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2200 { 2201 PetscFunctionBegin; 2202 PetscFunctionReturn(0); 2203 } 2204 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2207 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2208 { 2209 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2210 PetscErrorCode ierr; 2211 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2212 PetscScalar dx,*y,*xx,*w3_array; 2213 PetscScalar *vscale_array; 2214 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2215 Vec w1,w2,w3; 2216 void *fctx = coloring->fctx; 2217 PetscBool flg = PETSC_FALSE; 2218 2219 PetscFunctionBegin; 2220 if (!coloring->w1) { 2221 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2222 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2223 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2224 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2225 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2226 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2227 } 2228 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2229 2230 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2231 ierr = PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2232 if (flg) { 2233 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2234 } else { 2235 PetscBool assembled; 2236 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2237 if (assembled) { 2238 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2239 } 2240 } 2241 2242 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2243 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2244 2245 /* 2246 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2247 coloring->F for the coarser grids from the finest 2248 */ 2249 if (coloring->F) { 2250 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2251 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2252 if (m1 != m2) { 2253 coloring->F = 0; 2254 } 2255 } 2256 2257 if (coloring->F) { 2258 w1 = coloring->F; 2259 coloring->F = 0; 2260 } else { 2261 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2262 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2263 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2264 } 2265 2266 /* 2267 Compute all the scale factors and share with other processors 2268 */ 2269 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2270 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2271 for (k=0; k<coloring->ncolors; k++) { 2272 /* 2273 Loop over each column associated with color adding the 2274 perturbation to the vector w3. 2275 */ 2276 for (l=0; l<coloring->ncolumns[k]; l++) { 2277 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2278 dx = xx[col]; 2279 if (dx == 0.0) dx = 1.0; 2280 #if !defined(PETSC_USE_COMPLEX) 2281 if (dx < umin && dx >= 0.0) dx = umin; 2282 else if (dx < 0.0 && dx > -umin) dx = -umin; 2283 #else 2284 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2285 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2286 #endif 2287 dx *= epsilon; 2288 vscale_array[col] = 1.0/dx; 2289 } 2290 } 2291 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2292 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2293 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2294 2295 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2296 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2297 2298 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2299 else vscaleforrow = coloring->columnsforrow; 2300 2301 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2302 /* 2303 Loop over each color 2304 */ 2305 for (k=0; k<coloring->ncolors; k++) { 2306 coloring->currentcolor = k; 2307 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2308 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2309 /* 2310 Loop over each column associated with color adding the 2311 perturbation to the vector w3. 2312 */ 2313 for (l=0; l<coloring->ncolumns[k]; l++) { 2314 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2315 dx = xx[col]; 2316 if (dx == 0.0) dx = 1.0; 2317 #if !defined(PETSC_USE_COMPLEX) 2318 if (dx < umin && dx >= 0.0) dx = umin; 2319 else if (dx < 0.0 && dx > -umin) dx = -umin; 2320 #else 2321 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2322 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2323 #endif 2324 dx *= epsilon; 2325 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2326 w3_array[col] += dx; 2327 } 2328 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2329 2330 /* 2331 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2332 */ 2333 2334 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2335 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2336 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2337 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2338 2339 /* 2340 Loop over rows of vector, putting results into Jacobian matrix 2341 */ 2342 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2343 for (l=0; l<coloring->nrows[k]; l++) { 2344 row = coloring->rows[k][l]; 2345 col = coloring->columnsforrow[k][l]; 2346 y[row] *= vscale_array[vscaleforrow[k][l]]; 2347 srow = row + start; 2348 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2349 } 2350 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2351 } 2352 coloring->currentcolor = k; 2353 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2354 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2355 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2356 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2357 PetscFunctionReturn(0); 2358 } 2359 2360 /* 2361 Computes the number of nonzeros per row needed for preallocation when X and Y 2362 have different nonzero structure. 2363 */ 2364 #undef __FUNCT__ 2365 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2366 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt* nnz) 2367 { 2368 PetscInt i,m=Y->rmap->N; 2369 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2370 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2371 const PetscInt *xi = x->i,*yi = y->i; 2372 2373 PetscFunctionBegin; 2374 /* Set the number of nonzeros in the new matrix */ 2375 for(i=0; i<m; i++) { 2376 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2377 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2378 nnz[i] = 0; 2379 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2380 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2381 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2382 nnz[i]++; 2383 } 2384 for (; k<nzy; k++) nnz[i]++; 2385 } 2386 PetscFunctionReturn(0); 2387 } 2388 2389 #undef __FUNCT__ 2390 #define __FUNCT__ "MatAXPY_SeqAIJ" 2391 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2392 { 2393 PetscErrorCode ierr; 2394 PetscInt i; 2395 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2396 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2397 2398 PetscFunctionBegin; 2399 if (str == SAME_NONZERO_PATTERN) { 2400 PetscScalar alpha = a; 2401 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2402 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2403 if (y->xtoy && y->XtoY != X) { 2404 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2405 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2406 } 2407 if (!y->xtoy) { /* get xtoy */ 2408 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2409 y->XtoY = X; 2410 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2411 } 2412 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2413 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2414 } else { 2415 Mat B; 2416 PetscInt *nnz; 2417 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2418 ierr = MatCreate(((PetscObject)Y)->comm,&B);CHKERRQ(ierr); 2419 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2420 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2421 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2422 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2423 ierr = MatSeqAIJSetPreallocation(B,PETSC_NULL,nnz);CHKERRQ(ierr); 2424 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2425 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2426 ierr = PetscFree(nnz);CHKERRQ(ierr); 2427 } 2428 PetscFunctionReturn(0); 2429 } 2430 2431 #undef __FUNCT__ 2432 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2433 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2434 { 2435 PetscErrorCode ierr; 2436 2437 PetscFunctionBegin; 2438 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2439 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2440 PetscFunctionReturn(0); 2441 } 2442 2443 #undef __FUNCT__ 2444 #define __FUNCT__ "MatConjugate_SeqAIJ" 2445 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2446 { 2447 #if defined(PETSC_USE_COMPLEX) 2448 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2449 PetscInt i,nz; 2450 PetscScalar *a; 2451 2452 PetscFunctionBegin; 2453 nz = aij->nz; 2454 a = aij->a; 2455 for (i=0; i<nz; i++) { 2456 a[i] = PetscConj(a[i]); 2457 } 2458 #else 2459 PetscFunctionBegin; 2460 #endif 2461 PetscFunctionReturn(0); 2462 } 2463 2464 #undef __FUNCT__ 2465 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2466 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2467 { 2468 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2469 PetscErrorCode ierr; 2470 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2471 PetscReal atmp; 2472 PetscScalar *x; 2473 MatScalar *aa; 2474 2475 PetscFunctionBegin; 2476 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2477 aa = a->a; 2478 ai = a->i; 2479 aj = a->j; 2480 2481 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2482 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2483 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2484 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2485 for (i=0; i<m; i++) { 2486 ncols = ai[1] - ai[0]; ai++; 2487 x[i] = 0.0; 2488 for (j=0; j<ncols; j++){ 2489 atmp = PetscAbsScalar(*aa); 2490 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2491 aa++; aj++; 2492 } 2493 } 2494 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 #undef __FUNCT__ 2499 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2500 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2501 { 2502 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2503 PetscErrorCode ierr; 2504 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2505 PetscScalar *x; 2506 MatScalar *aa; 2507 2508 PetscFunctionBegin; 2509 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2510 aa = a->a; 2511 ai = a->i; 2512 aj = a->j; 2513 2514 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2515 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2516 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2517 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2518 for (i=0; i<m; i++) { 2519 ncols = ai[1] - ai[0]; ai++; 2520 if (ncols == A->cmap->n) { /* row is dense */ 2521 x[i] = *aa; if (idx) idx[i] = 0; 2522 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2523 x[i] = 0.0; 2524 if (idx) { 2525 idx[i] = 0; /* in case ncols is zero */ 2526 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2527 if (aj[j] > j) { 2528 idx[i] = j; 2529 break; 2530 } 2531 } 2532 } 2533 } 2534 for (j=0; j<ncols; j++){ 2535 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2536 aa++; aj++; 2537 } 2538 } 2539 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2540 PetscFunctionReturn(0); 2541 } 2542 2543 #undef __FUNCT__ 2544 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2545 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2546 { 2547 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2548 PetscErrorCode ierr; 2549 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2550 PetscReal atmp; 2551 PetscScalar *x; 2552 MatScalar *aa; 2553 2554 PetscFunctionBegin; 2555 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2556 aa = a->a; 2557 ai = a->i; 2558 aj = a->j; 2559 2560 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2561 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2562 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2563 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2564 for (i=0; i<m; i++) { 2565 ncols = ai[1] - ai[0]; ai++; 2566 if (ncols) { 2567 /* Get first nonzero */ 2568 for(j = 0; j < ncols; j++) { 2569 atmp = PetscAbsScalar(aa[j]); 2570 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2571 } 2572 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2573 } else { 2574 x[i] = 0.0; if (idx) idx[i] = 0; 2575 } 2576 for(j = 0; j < ncols; j++) { 2577 atmp = PetscAbsScalar(*aa); 2578 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2579 aa++; aj++; 2580 } 2581 } 2582 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2583 PetscFunctionReturn(0); 2584 } 2585 2586 #undef __FUNCT__ 2587 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2588 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2589 { 2590 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2591 PetscErrorCode ierr; 2592 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2593 PetscScalar *x; 2594 MatScalar *aa; 2595 2596 PetscFunctionBegin; 2597 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2598 aa = a->a; 2599 ai = a->i; 2600 aj = a->j; 2601 2602 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2603 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2604 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2605 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2606 for (i=0; i<m; i++) { 2607 ncols = ai[1] - ai[0]; ai++; 2608 if (ncols == A->cmap->n) { /* row is dense */ 2609 x[i] = *aa; if (idx) idx[i] = 0; 2610 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2611 x[i] = 0.0; 2612 if (idx) { /* find first implicit 0.0 in the row */ 2613 idx[i] = 0; /* in case ncols is zero */ 2614 for (j=0;j<ncols;j++) { 2615 if (aj[j] > j) { 2616 idx[i] = j; 2617 break; 2618 } 2619 } 2620 } 2621 } 2622 for (j=0; j<ncols; j++){ 2623 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2624 aa++; aj++; 2625 } 2626 } 2627 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2628 PetscFunctionReturn(0); 2629 } 2630 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2631 /* -------------------------------------------------------------------*/ 2632 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2633 MatGetRow_SeqAIJ, 2634 MatRestoreRow_SeqAIJ, 2635 MatMult_SeqAIJ, 2636 /* 4*/ MatMultAdd_SeqAIJ, 2637 MatMultTranspose_SeqAIJ, 2638 MatMultTransposeAdd_SeqAIJ, 2639 0, 2640 0, 2641 0, 2642 /*10*/ 0, 2643 MatLUFactor_SeqAIJ, 2644 0, 2645 MatSOR_SeqAIJ, 2646 MatTranspose_SeqAIJ, 2647 /*15*/ MatGetInfo_SeqAIJ, 2648 MatEqual_SeqAIJ, 2649 MatGetDiagonal_SeqAIJ, 2650 MatDiagonalScale_SeqAIJ, 2651 MatNorm_SeqAIJ, 2652 /*20*/ 0, 2653 MatAssemblyEnd_SeqAIJ, 2654 MatSetOption_SeqAIJ, 2655 MatZeroEntries_SeqAIJ, 2656 /*24*/ MatZeroRows_SeqAIJ, 2657 0, 2658 0, 2659 0, 2660 0, 2661 /*29*/ MatSetUpPreallocation_SeqAIJ, 2662 0, 2663 0, 2664 MatGetArray_SeqAIJ, 2665 MatRestoreArray_SeqAIJ, 2666 /*34*/ MatDuplicate_SeqAIJ, 2667 0, 2668 0, 2669 MatILUFactor_SeqAIJ, 2670 0, 2671 /*39*/ MatAXPY_SeqAIJ, 2672 MatGetSubMatrices_SeqAIJ, 2673 MatIncreaseOverlap_SeqAIJ, 2674 MatGetValues_SeqAIJ, 2675 MatCopy_SeqAIJ, 2676 /*44*/ MatGetRowMax_SeqAIJ, 2677 MatScale_SeqAIJ, 2678 0, 2679 MatDiagonalSet_SeqAIJ, 2680 MatZeroRowsColumns_SeqAIJ, 2681 /*49*/ MatSetBlockSize_SeqAIJ, 2682 MatGetRowIJ_SeqAIJ, 2683 MatRestoreRowIJ_SeqAIJ, 2684 MatGetColumnIJ_SeqAIJ, 2685 MatRestoreColumnIJ_SeqAIJ, 2686 /*54*/ MatFDColoringCreate_SeqAIJ, 2687 0, 2688 0, 2689 MatPermute_SeqAIJ, 2690 0, 2691 /*59*/ 0, 2692 MatDestroy_SeqAIJ, 2693 MatView_SeqAIJ, 2694 0, 2695 0, 2696 /*64*/ 0, 2697 0, 2698 0, 2699 0, 2700 0, 2701 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2702 MatGetRowMinAbs_SeqAIJ, 2703 0, 2704 MatSetColoring_SeqAIJ, 2705 #if defined(PETSC_HAVE_ADIC) 2706 MatSetValuesAdic_SeqAIJ, 2707 #else 2708 0, 2709 #endif 2710 /*74*/ MatSetValuesAdifor_SeqAIJ, 2711 MatFDColoringApply_AIJ, 2712 0, 2713 0, 2714 0, 2715 /*79*/ 0, 2716 0, 2717 0, 2718 0, 2719 MatLoad_SeqAIJ, 2720 /*84*/ MatIsSymmetric_SeqAIJ, 2721 MatIsHermitian_SeqAIJ, 2722 0, 2723 0, 2724 0, 2725 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2726 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2727 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2728 MatPtAP_Basic, 2729 MatPtAPSymbolic_SeqAIJ, 2730 /*94*/ MatPtAPNumeric_SeqAIJ, 2731 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2732 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2733 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2734 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2735 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2736 0, 2737 0, 2738 MatConjugate_SeqAIJ, 2739 0, 2740 /*104*/MatSetValuesRow_SeqAIJ, 2741 MatRealPart_SeqAIJ, 2742 MatImaginaryPart_SeqAIJ, 2743 0, 2744 0, 2745 /*109*/0, 2746 0, 2747 MatGetRowMin_SeqAIJ, 2748 0, 2749 MatMissingDiagonal_SeqAIJ, 2750 /*114*/0, 2751 0, 2752 0, 2753 0, 2754 0, 2755 /*119*/0, 2756 0, 2757 0, 2758 0, 2759 MatGetMultiProcBlock_SeqAIJ 2760 }; 2761 2762 EXTERN_C_BEGIN 2763 #undef __FUNCT__ 2764 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2765 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2766 { 2767 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2768 PetscInt i,nz,n; 2769 2770 PetscFunctionBegin; 2771 2772 nz = aij->maxnz; 2773 n = mat->rmap->n; 2774 for (i=0; i<nz; i++) { 2775 aij->j[i] = indices[i]; 2776 } 2777 aij->nz = nz; 2778 for (i=0; i<n; i++) { 2779 aij->ilen[i] = aij->imax[i]; 2780 } 2781 2782 PetscFunctionReturn(0); 2783 } 2784 EXTERN_C_END 2785 2786 #undef __FUNCT__ 2787 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2788 /*@ 2789 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2790 in the matrix. 2791 2792 Input Parameters: 2793 + mat - the SeqAIJ matrix 2794 - indices - the column indices 2795 2796 Level: advanced 2797 2798 Notes: 2799 This can be called if you have precomputed the nonzero structure of the 2800 matrix and want to provide it to the matrix object to improve the performance 2801 of the MatSetValues() operation. 2802 2803 You MUST have set the correct numbers of nonzeros per row in the call to 2804 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2805 2806 MUST be called before any calls to MatSetValues(); 2807 2808 The indices should start with zero, not one. 2809 2810 @*/ 2811 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2812 { 2813 PetscErrorCode ierr; 2814 2815 PetscFunctionBegin; 2816 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2817 PetscValidPointer(indices,2); 2818 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 2819 PetscFunctionReturn(0); 2820 } 2821 2822 /* ----------------------------------------------------------------------------------------*/ 2823 2824 EXTERN_C_BEGIN 2825 #undef __FUNCT__ 2826 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2827 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2828 { 2829 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2830 PetscErrorCode ierr; 2831 size_t nz = aij->i[mat->rmap->n]; 2832 2833 PetscFunctionBegin; 2834 if (aij->nonew != 1) { 2835 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2836 } 2837 2838 /* allocate space for values if not already there */ 2839 if (!aij->saved_values) { 2840 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2841 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2842 } 2843 2844 /* copy values over */ 2845 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2846 PetscFunctionReturn(0); 2847 } 2848 EXTERN_C_END 2849 2850 #undef __FUNCT__ 2851 #define __FUNCT__ "MatStoreValues" 2852 /*@ 2853 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2854 example, reuse of the linear part of a Jacobian, while recomputing the 2855 nonlinear portion. 2856 2857 Collect on Mat 2858 2859 Input Parameters: 2860 . mat - the matrix (currently only AIJ matrices support this option) 2861 2862 Level: advanced 2863 2864 Common Usage, with SNESSolve(): 2865 $ Create Jacobian matrix 2866 $ Set linear terms into matrix 2867 $ Apply boundary conditions to matrix, at this time matrix must have 2868 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2869 $ boundary conditions again will not change the nonzero structure 2870 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2871 $ ierr = MatStoreValues(mat); 2872 $ Call SNESSetJacobian() with matrix 2873 $ In your Jacobian routine 2874 $ ierr = MatRetrieveValues(mat); 2875 $ Set nonlinear terms in matrix 2876 2877 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2878 $ // build linear portion of Jacobian 2879 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2880 $ ierr = MatStoreValues(mat); 2881 $ loop over nonlinear iterations 2882 $ ierr = MatRetrieveValues(mat); 2883 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2884 $ // call MatAssemblyBegin/End() on matrix 2885 $ Solve linear system with Jacobian 2886 $ endloop 2887 2888 Notes: 2889 Matrix must already be assemblied before calling this routine 2890 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2891 calling this routine. 2892 2893 When this is called multiple times it overwrites the previous set of stored values 2894 and does not allocated additional space. 2895 2896 .seealso: MatRetrieveValues() 2897 2898 @*/ 2899 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2900 { 2901 PetscErrorCode ierr; 2902 2903 PetscFunctionBegin; 2904 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2905 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2906 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2907 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 2908 PetscFunctionReturn(0); 2909 } 2910 2911 EXTERN_C_BEGIN 2912 #undef __FUNCT__ 2913 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2914 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2915 { 2916 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2917 PetscErrorCode ierr; 2918 PetscInt nz = aij->i[mat->rmap->n]; 2919 2920 PetscFunctionBegin; 2921 if (aij->nonew != 1) { 2922 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2923 } 2924 if (!aij->saved_values) { 2925 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2926 } 2927 /* copy values over */ 2928 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2929 PetscFunctionReturn(0); 2930 } 2931 EXTERN_C_END 2932 2933 #undef __FUNCT__ 2934 #define __FUNCT__ "MatRetrieveValues" 2935 /*@ 2936 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2937 example, reuse of the linear part of a Jacobian, while recomputing the 2938 nonlinear portion. 2939 2940 Collect on Mat 2941 2942 Input Parameters: 2943 . mat - the matrix (currently on AIJ matrices support this option) 2944 2945 Level: advanced 2946 2947 .seealso: MatStoreValues() 2948 2949 @*/ 2950 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2951 { 2952 PetscErrorCode ierr; 2953 2954 PetscFunctionBegin; 2955 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2956 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2957 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2958 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 2959 PetscFunctionReturn(0); 2960 } 2961 2962 2963 /* --------------------------------------------------------------------------------*/ 2964 #undef __FUNCT__ 2965 #define __FUNCT__ "MatCreateSeqAIJ" 2966 /*@C 2967 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2968 (the default parallel PETSc format). For good matrix assembly performance 2969 the user should preallocate the matrix storage by setting the parameter nz 2970 (or the array nnz). By setting these parameters accurately, performance 2971 during matrix assembly can be increased by more than a factor of 50. 2972 2973 Collective on MPI_Comm 2974 2975 Input Parameters: 2976 + comm - MPI communicator, set to PETSC_COMM_SELF 2977 . m - number of rows 2978 . n - number of columns 2979 . nz - number of nonzeros per row (same for all rows) 2980 - nnz - array containing the number of nonzeros in the various rows 2981 (possibly different for each row) or PETSC_NULL 2982 2983 Output Parameter: 2984 . A - the matrix 2985 2986 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2987 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2988 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2989 2990 Notes: 2991 If nnz is given then nz is ignored 2992 2993 The AIJ format (also called the Yale sparse matrix format or 2994 compressed row storage), is fully compatible with standard Fortran 77 2995 storage. That is, the stored row and column indices can begin at 2996 either one (as in Fortran) or zero. See the users' manual for details. 2997 2998 Specify the preallocated storage with either nz or nnz (not both). 2999 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3000 allocation. For large problems you MUST preallocate memory or you 3001 will get TERRIBLE performance, see the users' manual chapter on matrices. 3002 3003 By default, this format uses inodes (identical nodes) when possible, to 3004 improve numerical efficiency of matrix-vector products and solves. We 3005 search for consecutive rows with the same nonzero structure, thereby 3006 reusing matrix information to achieve increased efficiency. 3007 3008 Options Database Keys: 3009 + -mat_no_inode - Do not use inodes 3010 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3011 3012 Level: intermediate 3013 3014 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3015 3016 @*/ 3017 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3018 { 3019 PetscErrorCode ierr; 3020 3021 PetscFunctionBegin; 3022 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3023 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3024 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3025 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3026 PetscFunctionReturn(0); 3027 } 3028 3029 #undef __FUNCT__ 3030 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3031 /*@C 3032 MatSeqAIJSetPreallocation - For good matrix assembly performance 3033 the user should preallocate the matrix storage by setting the parameter nz 3034 (or the array nnz). By setting these parameters accurately, performance 3035 during matrix assembly can be increased by more than a factor of 50. 3036 3037 Collective on MPI_Comm 3038 3039 Input Parameters: 3040 + B - The matrix-free 3041 . nz - number of nonzeros per row (same for all rows) 3042 - nnz - array containing the number of nonzeros in the various rows 3043 (possibly different for each row) or PETSC_NULL 3044 3045 Notes: 3046 If nnz is given then nz is ignored 3047 3048 The AIJ format (also called the Yale sparse matrix format or 3049 compressed row storage), is fully compatible with standard Fortran 77 3050 storage. That is, the stored row and column indices can begin at 3051 either one (as in Fortran) or zero. See the users' manual for details. 3052 3053 Specify the preallocated storage with either nz or nnz (not both). 3054 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3055 allocation. For large problems you MUST preallocate memory or you 3056 will get TERRIBLE performance, see the users' manual chapter on matrices. 3057 3058 You can call MatGetInfo() to get information on how effective the preallocation was; 3059 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3060 You can also run with the option -info and look for messages with the string 3061 malloc in them to see if additional memory allocation was needed. 3062 3063 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3064 entries or columns indices 3065 3066 By default, this format uses inodes (identical nodes) when possible, to 3067 improve numerical efficiency of matrix-vector products and solves. We 3068 search for consecutive rows with the same nonzero structure, thereby 3069 reusing matrix information to achieve increased efficiency. 3070 3071 Options Database Keys: 3072 + -mat_no_inode - Do not use inodes 3073 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3074 - -mat_aij_oneindex - Internally use indexing starting at 1 3075 rather than 0. Note that when calling MatSetValues(), 3076 the user still MUST index entries starting at 0! 3077 3078 Level: intermediate 3079 3080 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3081 3082 @*/ 3083 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3084 { 3085 PetscErrorCode ierr; 3086 3087 PetscFunctionBegin; 3088 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3089 PetscFunctionReturn(0); 3090 } 3091 3092 EXTERN_C_BEGIN 3093 #undef __FUNCT__ 3094 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3095 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3096 { 3097 Mat_SeqAIJ *b; 3098 PetscBool skipallocation = PETSC_FALSE; 3099 PetscErrorCode ierr; 3100 PetscInt i; 3101 3102 PetscFunctionBegin; 3103 3104 if (nz == MAT_SKIP_ALLOCATION) { 3105 skipallocation = PETSC_TRUE; 3106 nz = 0; 3107 } 3108 3109 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3110 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3111 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3112 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3113 3114 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3115 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3116 if (nnz) { 3117 for (i=0; i<B->rmap->n; i++) { 3118 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]); 3119 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); 3120 } 3121 } 3122 3123 B->preallocated = PETSC_TRUE; 3124 b = (Mat_SeqAIJ*)B->data; 3125 3126 if (!skipallocation) { 3127 if (!b->imax) { 3128 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3129 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3130 } 3131 if (!nnz) { 3132 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3133 else if (nz < 0) nz = 1; 3134 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3135 nz = nz*B->rmap->n; 3136 } else { 3137 nz = 0; 3138 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3139 } 3140 /* b->ilen will count nonzeros in each row so far. */ 3141 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3142 3143 /* allocate the matrix space */ 3144 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3145 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3146 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3147 b->i[0] = 0; 3148 for (i=1; i<B->rmap->n+1; i++) { 3149 b->i[i] = b->i[i-1] + b->imax[i-1]; 3150 } 3151 b->singlemalloc = PETSC_TRUE; 3152 b->free_a = PETSC_TRUE; 3153 b->free_ij = PETSC_TRUE; 3154 } else { 3155 b->free_a = PETSC_FALSE; 3156 b->free_ij = PETSC_FALSE; 3157 } 3158 3159 b->nz = 0; 3160 b->maxnz = nz; 3161 B->info.nz_unneeded = (double)b->maxnz; 3162 PetscFunctionReturn(0); 3163 } 3164 EXTERN_C_END 3165 3166 #undef __FUNCT__ 3167 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3168 /*@ 3169 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3170 3171 Input Parameters: 3172 + B - the matrix 3173 . i - the indices into j for the start of each row (starts with zero) 3174 . j - the column indices for each row (starts with zero) these must be sorted for each row 3175 - v - optional values in the matrix 3176 3177 Level: developer 3178 3179 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3180 3181 .keywords: matrix, aij, compressed row, sparse, sequential 3182 3183 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3184 @*/ 3185 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3186 { 3187 PetscErrorCode ierr; 3188 3189 PetscFunctionBegin; 3190 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3191 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3192 PetscFunctionReturn(0); 3193 } 3194 3195 EXTERN_C_BEGIN 3196 #undef __FUNCT__ 3197 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3198 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3199 { 3200 PetscInt i; 3201 PetscInt m,n; 3202 PetscInt nz; 3203 PetscInt *nnz, nz_max = 0; 3204 PetscScalar *values; 3205 PetscErrorCode ierr; 3206 3207 PetscFunctionBegin; 3208 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3209 3210 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3211 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3212 for(i = 0; i < m; i++) { 3213 nz = Ii[i+1]- Ii[i]; 3214 nz_max = PetscMax(nz_max, nz); 3215 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3216 nnz[i] = nz; 3217 } 3218 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3219 ierr = PetscFree(nnz);CHKERRQ(ierr); 3220 3221 if (v) { 3222 values = (PetscScalar*) v; 3223 } else { 3224 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3225 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3226 } 3227 3228 for(i = 0; i < m; i++) { 3229 nz = Ii[i+1] - Ii[i]; 3230 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3231 } 3232 3233 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3234 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3235 3236 if (!v) { 3237 ierr = PetscFree(values);CHKERRQ(ierr); 3238 } 3239 PetscFunctionReturn(0); 3240 } 3241 EXTERN_C_END 3242 3243 #include "../src/mat/impls/dense/seq/dense.h" 3244 #include "private/petscaxpy.h" 3245 3246 #undef __FUNCT__ 3247 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3248 /* 3249 Computes (B'*A')' since computing B*A directly is untenable 3250 3251 n p p 3252 ( ) ( ) ( ) 3253 m ( A ) * n ( B ) = m ( C ) 3254 ( ) ( ) ( ) 3255 3256 */ 3257 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3258 { 3259 PetscErrorCode ierr; 3260 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3261 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3262 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3263 PetscInt i,n,m,q,p; 3264 const PetscInt *ii,*idx; 3265 const PetscScalar *b,*a,*a_q; 3266 PetscScalar *c,*c_q; 3267 3268 PetscFunctionBegin; 3269 m = A->rmap->n; 3270 n = A->cmap->n; 3271 p = B->cmap->n; 3272 a = sub_a->v; 3273 b = sub_b->a; 3274 c = sub_c->v; 3275 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3276 3277 ii = sub_b->i; 3278 idx = sub_b->j; 3279 for (i=0; i<n; i++) { 3280 q = ii[i+1] - ii[i]; 3281 while (q-->0) { 3282 c_q = c + m*(*idx); 3283 a_q = a + m*i; 3284 PetscAXPY(c_q,*b,a_q,m); 3285 idx++; 3286 b++; 3287 } 3288 } 3289 PetscFunctionReturn(0); 3290 } 3291 3292 #undef __FUNCT__ 3293 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3294 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3295 { 3296 PetscErrorCode ierr; 3297 PetscInt m=A->rmap->n,n=B->cmap->n; 3298 Mat Cmat; 3299 3300 PetscFunctionBegin; 3301 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); 3302 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3303 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3304 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3305 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3306 Cmat->assembled = PETSC_TRUE; 3307 *C = Cmat; 3308 PetscFunctionReturn(0); 3309 } 3310 3311 /* ----------------------------------------------------------------*/ 3312 #undef __FUNCT__ 3313 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3314 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3315 { 3316 PetscErrorCode ierr; 3317 3318 PetscFunctionBegin; 3319 if (scall == MAT_INITIAL_MATRIX){ 3320 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3321 } 3322 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3323 PetscFunctionReturn(0); 3324 } 3325 3326 3327 /*MC 3328 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3329 based on compressed sparse row format. 3330 3331 Options Database Keys: 3332 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3333 3334 Level: beginner 3335 3336 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3337 M*/ 3338 3339 EXTERN_C_BEGIN 3340 #if defined(PETSC_HAVE_PASTIX) 3341 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3342 #endif 3343 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3344 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3345 #endif 3346 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3347 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3348 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3349 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3350 #if defined(PETSC_HAVE_MUMPS) 3351 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3352 #endif 3353 #if defined(PETSC_HAVE_SUPERLU) 3354 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3355 #endif 3356 #if defined(PETSC_HAVE_SUPERLU_DIST) 3357 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3358 #endif 3359 #if defined(PETSC_HAVE_SPOOLES) 3360 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3361 #endif 3362 #if defined(PETSC_HAVE_UMFPACK) 3363 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3364 #endif 3365 #if defined(PETSC_HAVE_CHOLMOD) 3366 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3367 #endif 3368 #if defined(PETSC_HAVE_LUSOL) 3369 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3370 #endif 3371 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3372 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3373 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3374 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3375 #endif 3376 EXTERN_C_END 3377 3378 3379 EXTERN_C_BEGIN 3380 #undef __FUNCT__ 3381 #define __FUNCT__ "MatCreate_SeqAIJ" 3382 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3383 { 3384 Mat_SeqAIJ *b; 3385 PetscErrorCode ierr; 3386 PetscMPIInt size; 3387 3388 PetscFunctionBegin; 3389 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3390 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3391 3392 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3393 B->data = (void*)b; 3394 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3395 B->mapping = 0; 3396 b->row = 0; 3397 b->col = 0; 3398 b->icol = 0; 3399 b->reallocs = 0; 3400 b->ignorezeroentries = PETSC_FALSE; 3401 b->roworiented = PETSC_TRUE; 3402 b->nonew = 0; 3403 b->diag = 0; 3404 b->solve_work = 0; 3405 B->spptr = 0; 3406 b->saved_values = 0; 3407 b->idiag = 0; 3408 b->mdiag = 0; 3409 b->ssor_work = 0; 3410 b->omega = 1.0; 3411 b->fshift = 0.0; 3412 b->idiagvalid = PETSC_FALSE; 3413 b->keepnonzeropattern = PETSC_FALSE; 3414 b->xtoy = 0; 3415 b->XtoY = 0; 3416 b->compressedrow.use = PETSC_FALSE; 3417 b->compressedrow.nrows = B->rmap->n; 3418 b->compressedrow.i = PETSC_NULL; 3419 b->compressedrow.rindex = PETSC_NULL; 3420 b->compressedrow.checked = PETSC_FALSE; 3421 B->same_nonzero = PETSC_FALSE; 3422 3423 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3424 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3425 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3426 "MatGetFactor_seqaij_matlab", 3427 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3428 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3429 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3430 #endif 3431 #if defined(PETSC_HAVE_PASTIX) 3432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3433 "MatGetFactor_seqaij_pastix", 3434 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3435 #endif 3436 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3437 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3438 "MatGetFactor_seqaij_essl", 3439 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3440 #endif 3441 #if defined(PETSC_HAVE_SUPERLU) 3442 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3443 "MatGetFactor_seqaij_superlu", 3444 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3445 #endif 3446 #if defined(PETSC_HAVE_SUPERLU_DIST) 3447 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3448 "MatGetFactor_seqaij_superlu_dist", 3449 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3450 #endif 3451 #if defined(PETSC_HAVE_SPOOLES) 3452 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3453 "MatGetFactor_seqaij_spooles", 3454 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3455 #endif 3456 #if defined(PETSC_HAVE_MUMPS) 3457 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3458 "MatGetFactor_aij_mumps", 3459 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3460 #endif 3461 #if defined(PETSC_HAVE_UMFPACK) 3462 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3463 "MatGetFactor_seqaij_umfpack", 3464 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3465 #endif 3466 #if defined(PETSC_HAVE_CHOLMOD) 3467 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3468 "MatGetFactor_seqaij_cholmod", 3469 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3470 #endif 3471 #if defined(PETSC_HAVE_LUSOL) 3472 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3473 "MatGetFactor_seqaij_lusol", 3474 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3475 #endif 3476 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3477 "MatGetFactor_seqaij_petsc", 3478 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3479 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3480 "MatGetFactorAvailable_seqaij_petsc", 3481 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3482 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3483 "MatGetFactor_seqaij_bas", 3484 MatGetFactor_seqaij_bas); 3485 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3486 "MatSeqAIJSetColumnIndices_SeqAIJ", 3487 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3488 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3489 "MatStoreValues_SeqAIJ", 3490 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3491 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3492 "MatRetrieveValues_SeqAIJ", 3493 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3494 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3495 "MatConvert_SeqAIJ_SeqSBAIJ", 3496 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3497 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3498 "MatConvert_SeqAIJ_SeqBAIJ", 3499 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3500 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C", 3501 "MatConvert_SeqAIJ_SeqAIJPERM", 3502 MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3503 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C", 3504 "MatConvert_SeqAIJ_SeqAIJCRL", 3505 MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3506 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3507 "MatIsTranspose_SeqAIJ", 3508 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3509 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3510 "MatIsHermitianTranspose_SeqAIJ", 3511 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3512 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3513 "MatSeqAIJSetPreallocation_SeqAIJ", 3514 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3515 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3516 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3517 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3518 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3519 "MatReorderForNonzeroDiagonal_SeqAIJ", 3520 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3521 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3522 "MatMatMult_SeqDense_SeqAIJ", 3523 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3524 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3525 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3526 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3527 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3528 "MatMatMultNumeric_SeqDense_SeqAIJ", 3529 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3530 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3531 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3532 PetscFunctionReturn(0); 3533 } 3534 EXTERN_C_END 3535 3536 #undef __FUNCT__ 3537 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3538 /* 3539 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3540 */ 3541 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3542 { 3543 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3544 PetscErrorCode ierr; 3545 PetscInt i,m = A->rmap->n; 3546 3547 PetscFunctionBegin; 3548 c = (Mat_SeqAIJ*)C->data; 3549 3550 C->factortype = A->factortype; 3551 c->row = 0; 3552 c->col = 0; 3553 c->icol = 0; 3554 c->reallocs = 0; 3555 3556 C->assembled = PETSC_TRUE; 3557 3558 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3559 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3560 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3561 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3562 3563 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3564 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3565 for (i=0; i<m; i++) { 3566 c->imax[i] = a->imax[i]; 3567 c->ilen[i] = a->ilen[i]; 3568 } 3569 3570 /* allocate the matrix space */ 3571 if (mallocmatspace){ 3572 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3573 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3574 c->singlemalloc = PETSC_TRUE; 3575 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3576 if (m > 0) { 3577 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3578 if (cpvalues == MAT_COPY_VALUES) { 3579 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3580 } else { 3581 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3582 } 3583 } 3584 } 3585 3586 c->ignorezeroentries = a->ignorezeroentries; 3587 c->roworiented = a->roworiented; 3588 c->nonew = a->nonew; 3589 if (a->diag) { 3590 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3591 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3592 for (i=0; i<m; i++) { 3593 c->diag[i] = a->diag[i]; 3594 } 3595 } else c->diag = 0; 3596 c->solve_work = 0; 3597 c->saved_values = 0; 3598 c->idiag = 0; 3599 c->ssor_work = 0; 3600 c->keepnonzeropattern = a->keepnonzeropattern; 3601 c->free_a = PETSC_TRUE; 3602 c->free_ij = PETSC_TRUE; 3603 c->xtoy = 0; 3604 c->XtoY = 0; 3605 3606 c->nz = a->nz; 3607 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3608 C->preallocated = PETSC_TRUE; 3609 3610 c->compressedrow.use = a->compressedrow.use; 3611 c->compressedrow.nrows = a->compressedrow.nrows; 3612 c->compressedrow.checked = a->compressedrow.checked; 3613 if (a->compressedrow.checked && a->compressedrow.use){ 3614 i = a->compressedrow.nrows; 3615 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3616 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3617 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3618 } else { 3619 c->compressedrow.use = PETSC_FALSE; 3620 c->compressedrow.i = PETSC_NULL; 3621 c->compressedrow.rindex = PETSC_NULL; 3622 } 3623 C->same_nonzero = A->same_nonzero; 3624 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3625 3626 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3627 PetscFunctionReturn(0); 3628 } 3629 3630 #undef __FUNCT__ 3631 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3632 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3633 { 3634 PetscErrorCode ierr; 3635 3636 PetscFunctionBegin; 3637 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3638 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3639 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3640 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3641 PetscFunctionReturn(0); 3642 } 3643 3644 #undef __FUNCT__ 3645 #define __FUNCT__ "MatLoad_SeqAIJ" 3646 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3647 { 3648 Mat_SeqAIJ *a; 3649 PetscErrorCode ierr; 3650 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3651 int fd; 3652 PetscMPIInt size; 3653 MPI_Comm comm; 3654 3655 PetscFunctionBegin; 3656 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3657 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3658 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3659 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3660 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3661 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3662 M = header[1]; N = header[2]; nz = header[3]; 3663 3664 if (nz < 0) { 3665 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3666 } 3667 3668 /* read in row lengths */ 3669 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3670 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3671 3672 /* check if sum of rowlengths is same as nz */ 3673 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3674 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); 3675 3676 /* set global size if not set already*/ 3677 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3678 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3679 } else { 3680 /* if sizes and type are already set, check if the vector global sizes are correct */ 3681 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3682 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); 3683 } 3684 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3685 a = (Mat_SeqAIJ*)newMat->data; 3686 3687 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3688 3689 /* read in nonzero values */ 3690 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3691 3692 /* set matrix "i" values */ 3693 a->i[0] = 0; 3694 for (i=1; i<= M; i++) { 3695 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3696 a->ilen[i-1] = rowlengths[i-1]; 3697 } 3698 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3699 3700 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3701 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3702 PetscFunctionReturn(0); 3703 } 3704 3705 #undef __FUNCT__ 3706 #define __FUNCT__ "MatEqual_SeqAIJ" 3707 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 3708 { 3709 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3710 PetscErrorCode ierr; 3711 #if defined(PETSC_USE_COMPLEX) 3712 PetscInt k; 3713 #endif 3714 3715 PetscFunctionBegin; 3716 /* If the matrix dimensions are not equal,or no of nonzeros */ 3717 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3718 *flg = PETSC_FALSE; 3719 PetscFunctionReturn(0); 3720 } 3721 3722 /* if the a->i are the same */ 3723 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3724 if (!*flg) PetscFunctionReturn(0); 3725 3726 /* if a->j are the same */ 3727 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3728 if (!*flg) PetscFunctionReturn(0); 3729 3730 /* if a->a are the same */ 3731 #if defined(PETSC_USE_COMPLEX) 3732 for (k=0; k<a->nz; k++){ 3733 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3734 *flg = PETSC_FALSE; 3735 PetscFunctionReturn(0); 3736 } 3737 } 3738 #else 3739 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3740 #endif 3741 PetscFunctionReturn(0); 3742 } 3743 3744 #undef __FUNCT__ 3745 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3746 /*@ 3747 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3748 provided by the user. 3749 3750 Collective on MPI_Comm 3751 3752 Input Parameters: 3753 + comm - must be an MPI communicator of size 1 3754 . m - number of rows 3755 . n - number of columns 3756 . i - row indices 3757 . j - column indices 3758 - a - matrix values 3759 3760 Output Parameter: 3761 . mat - the matrix 3762 3763 Level: intermediate 3764 3765 Notes: 3766 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3767 once the matrix is destroyed 3768 3769 You cannot set new nonzero locations into this matrix, that will generate an error. 3770 3771 The i and j indices are 0 based 3772 3773 The format which is used for the sparse matrix input, is equivalent to a 3774 row-major ordering.. i.e for the following matrix, the input data expected is 3775 as shown: 3776 3777 1 0 0 3778 2 0 3 3779 4 5 6 3780 3781 i = {0,1,3,6} [size = nrow+1 = 3+1] 3782 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3783 v = {1,2,3,4,5,6} [size = nz = 6] 3784 3785 3786 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3787 3788 @*/ 3789 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3790 { 3791 PetscErrorCode ierr; 3792 PetscInt ii; 3793 Mat_SeqAIJ *aij; 3794 #if defined(PETSC_USE_DEBUG) 3795 PetscInt jj; 3796 #endif 3797 3798 PetscFunctionBegin; 3799 if (i[0]) { 3800 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3801 } 3802 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3803 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3804 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3805 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3806 aij = (Mat_SeqAIJ*)(*mat)->data; 3807 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3808 3809 aij->i = i; 3810 aij->j = j; 3811 aij->a = a; 3812 aij->singlemalloc = PETSC_FALSE; 3813 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3814 aij->free_a = PETSC_FALSE; 3815 aij->free_ij = PETSC_FALSE; 3816 3817 for (ii=0; ii<m; ii++) { 3818 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3819 #if defined(PETSC_USE_DEBUG) 3820 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]); 3821 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3822 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); 3823 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); 3824 } 3825 #endif 3826 } 3827 #if defined(PETSC_USE_DEBUG) 3828 for (ii=0; ii<aij->i[m]; ii++) { 3829 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3830 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]); 3831 } 3832 #endif 3833 3834 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3835 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3836 PetscFunctionReturn(0); 3837 } 3838 3839 #undef __FUNCT__ 3840 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3841 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3842 { 3843 PetscErrorCode ierr; 3844 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3845 3846 PetscFunctionBegin; 3847 if (coloring->ctype == IS_COLORING_GLOBAL) { 3848 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3849 a->coloring = coloring; 3850 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3851 PetscInt i,*larray; 3852 ISColoring ocoloring; 3853 ISColoringValue *colors; 3854 3855 /* set coloring for diagonal portion */ 3856 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3857 for (i=0; i<A->cmap->n; i++) { 3858 larray[i] = i; 3859 } 3860 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3861 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3862 for (i=0; i<A->cmap->n; i++) { 3863 colors[i] = coloring->colors[larray[i]]; 3864 } 3865 ierr = PetscFree(larray);CHKERRQ(ierr); 3866 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3867 a->coloring = ocoloring; 3868 } 3869 PetscFunctionReturn(0); 3870 } 3871 3872 #if defined(PETSC_HAVE_ADIC) 3873 EXTERN_C_BEGIN 3874 #include "adic/ad_utils.h" 3875 EXTERN_C_END 3876 3877 #undef __FUNCT__ 3878 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3879 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3880 { 3881 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3882 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3883 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3884 ISColoringValue *color; 3885 3886 PetscFunctionBegin; 3887 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3888 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3889 color = a->coloring->colors; 3890 /* loop over rows */ 3891 for (i=0; i<m; i++) { 3892 nz = ii[i+1] - ii[i]; 3893 /* loop over columns putting computed value into matrix */ 3894 for (j=0; j<nz; j++) { 3895 *v++ = values[color[*jj++]]; 3896 } 3897 values += nlen; /* jump to next row of derivatives */ 3898 } 3899 PetscFunctionReturn(0); 3900 } 3901 #endif 3902 3903 #undef __FUNCT__ 3904 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3905 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3906 { 3907 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3908 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3909 MatScalar *v = a->a; 3910 PetscScalar *values = (PetscScalar *)advalues; 3911 ISColoringValue *color; 3912 3913 PetscFunctionBegin; 3914 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3915 color = a->coloring->colors; 3916 /* loop over rows */ 3917 for (i=0; i<m; i++) { 3918 nz = ii[i+1] - ii[i]; 3919 /* loop over columns putting computed value into matrix */ 3920 for (j=0; j<nz; j++) { 3921 *v++ = values[color[*jj++]]; 3922 } 3923 values += nl; /* jump to next row of derivatives */ 3924 } 3925 PetscFunctionReturn(0); 3926 } 3927 3928 /* 3929 Special version for direct calls from Fortran 3930 */ 3931 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3932 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3933 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3934 #define matsetvaluesseqaij_ matsetvaluesseqaij 3935 #endif 3936 3937 /* Change these macros so can be used in void function */ 3938 #undef CHKERRQ 3939 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3940 #undef SETERRQ2 3941 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3942 3943 EXTERN_C_BEGIN 3944 #undef __FUNCT__ 3945 #define __FUNCT__ "matsetvaluesseqaij_" 3946 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3947 { 3948 Mat A = *AA; 3949 PetscInt m = *mm, n = *nn; 3950 InsertMode is = *isis; 3951 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3952 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3953 PetscInt *imax,*ai,*ailen; 3954 PetscErrorCode ierr; 3955 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3956 MatScalar *ap,value,*aa; 3957 PetscBool ignorezeroentries = a->ignorezeroentries; 3958 PetscBool roworiented = a->roworiented; 3959 3960 PetscFunctionBegin; 3961 ierr = MatPreallocated(A);CHKERRQ(ierr); 3962 imax = a->imax; 3963 ai = a->i; 3964 ailen = a->ilen; 3965 aj = a->j; 3966 aa = a->a; 3967 3968 for (k=0; k<m; k++) { /* loop over added rows */ 3969 row = im[k]; 3970 if (row < 0) continue; 3971 #if defined(PETSC_USE_DEBUG) 3972 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3973 #endif 3974 rp = aj + ai[row]; ap = aa + ai[row]; 3975 rmax = imax[row]; nrow = ailen[row]; 3976 low = 0; 3977 high = nrow; 3978 for (l=0; l<n; l++) { /* loop over added columns */ 3979 if (in[l] < 0) continue; 3980 #if defined(PETSC_USE_DEBUG) 3981 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3982 #endif 3983 col = in[l]; 3984 if (roworiented) { 3985 value = v[l + k*n]; 3986 } else { 3987 value = v[k + l*m]; 3988 } 3989 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3990 3991 if (col <= lastcol) low = 0; else high = nrow; 3992 lastcol = col; 3993 while (high-low > 5) { 3994 t = (low+high)/2; 3995 if (rp[t] > col) high = t; 3996 else low = t; 3997 } 3998 for (i=low; i<high; i++) { 3999 if (rp[i] > col) break; 4000 if (rp[i] == col) { 4001 if (is == ADD_VALUES) ap[i] += value; 4002 else ap[i] = value; 4003 goto noinsert; 4004 } 4005 } 4006 if (value == 0.0 && ignorezeroentries) goto noinsert; 4007 if (nonew == 1) goto noinsert; 4008 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4009 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4010 N = nrow++ - 1; a->nz++; high++; 4011 /* shift up all the later entries in this row */ 4012 for (ii=N; ii>=i; ii--) { 4013 rp[ii+1] = rp[ii]; 4014 ap[ii+1] = ap[ii]; 4015 } 4016 rp[i] = col; 4017 ap[i] = value; 4018 noinsert:; 4019 low = i + 1; 4020 } 4021 ailen[row] = nrow; 4022 } 4023 A->same_nonzero = PETSC_FALSE; 4024 PetscFunctionReturnVoid(); 4025 } 4026 EXTERN_C_END 4027