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