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