1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the AIJ (compressed row) 5 matrix storage format. 6 */ 7 8 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 9 #include "src/inline/spops.h" 10 #include "src/inline/dot.h" 11 #include "petscbt.h" 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 15 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 16 { 17 PetscErrorCode ierr; 18 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 19 PetscInt i,*diag, m = Y->rmap.n; 20 PetscScalar *v,*aa = aij->a; 21 PetscTruth missing; 22 23 PetscFunctionBegin; 24 if (Y->assembled) { 25 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 26 if (!missing) { 27 diag = aij->diag; 28 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 29 if (is == INSERT_VALUES) { 30 for (i=0; i<m; i++) { 31 aa[diag[i]] = v[i]; 32 } 33 } else { 34 for (i=0; i<m; i++) { 35 aa[diag[i]] += v[i]; 36 } 37 } 38 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 39 PetscFunctionReturn(0); 40 } 41 } 42 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 48 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 49 { 50 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 51 PetscErrorCode ierr; 52 PetscInt i,ishift; 53 54 PetscFunctionBegin; 55 *m = A->rmap.n; 56 if (!ia) PetscFunctionReturn(0); 57 ishift = 0; 58 if (symmetric && !A->structurally_symmetric) { 59 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap.n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 60 } else if (oshift == 1) { 61 PetscInt nz = a->i[A->rmap.n]; 62 /* malloc space and add 1 to i and j indices */ 63 ierr = PetscMalloc((A->rmap.n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 64 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 65 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 66 for (i=0; i<A->rmap.n+1; i++) (*ia)[i] = a->i[i] + 1; 67 } else { 68 *ia = a->i; *ja = a->j; 69 } 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 75 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 76 { 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 if (!ia) PetscFunctionReturn(0); 81 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 82 ierr = PetscFree(*ia);CHKERRQ(ierr); 83 ierr = PetscFree(*ja);CHKERRQ(ierr); 84 } 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 90 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 91 { 92 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 93 PetscErrorCode ierr; 94 PetscInt i,*collengths,*cia,*cja,n = A->cmap.n,m = A->rmap.n; 95 PetscInt nz = a->i[m],row,*jj,mr,col; 96 97 PetscFunctionBegin; 98 *nn = n; 99 if (!ia) PetscFunctionReturn(0); 100 if (symmetric) { 101 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap.n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 102 } else { 103 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 104 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 105 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 106 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 107 jj = a->j; 108 for (i=0; i<nz; i++) { 109 collengths[jj[i]]++; 110 } 111 cia[0] = oshift; 112 for (i=0; i<n; i++) { 113 cia[i+1] = cia[i] + collengths[i]; 114 } 115 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 116 jj = a->j; 117 for (row=0; row<m; row++) { 118 mr = a->i[row+1] - a->i[row]; 119 for (i=0; i<mr; i++) { 120 col = *jj++; 121 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 122 } 123 } 124 ierr = PetscFree(collengths);CHKERRQ(ierr); 125 *ia = cia; *ja = cja; 126 } 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 132 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 133 { 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 if (!ia) PetscFunctionReturn(0); 138 139 ierr = PetscFree(*ia);CHKERRQ(ierr); 140 ierr = PetscFree(*ja);CHKERRQ(ierr); 141 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 147 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 148 { 149 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 150 PetscInt *ai = a->i; 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #define CHUNKSIZE 15 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "MatSetValues_SeqAIJ" 162 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 163 { 164 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 165 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 166 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 167 PetscErrorCode ierr; 168 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 169 PetscScalar *ap,value,*aa = a->a; 170 PetscTruth ignorezeroentries = a->ignorezeroentries; 171 PetscTruth roworiented = a->roworiented; 172 173 PetscFunctionBegin; 174 for (k=0; k<m; k++) { /* loop over added rows */ 175 row = im[k]; 176 if (row < 0) continue; 177 #if defined(PETSC_USE_DEBUG) 178 if (row >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.n-1); 179 #endif 180 rp = aj + ai[row]; ap = aa + ai[row]; 181 rmax = imax[row]; nrow = ailen[row]; 182 low = 0; 183 high = nrow; 184 for (l=0; l<n; l++) { /* loop over added columns */ 185 if (in[l] < 0) continue; 186 #if defined(PETSC_USE_DEBUG) 187 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 188 #endif 189 col = in[l]; 190 if (roworiented) { 191 value = v[l + k*n]; 192 } else { 193 value = v[k + l*m]; 194 } 195 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 196 197 if (col <= lastcol) low = 0; else high = nrow; 198 lastcol = col; 199 while (high-low > 5) { 200 t = (low+high)/2; 201 if (rp[t] > col) high = t; 202 else low = t; 203 } 204 for (i=low; i<high; i++) { 205 if (rp[i] > col) break; 206 if (rp[i] == col) { 207 if (is == ADD_VALUES) ap[i] += value; 208 else ap[i] = value; 209 goto noinsert; 210 } 211 } 212 if (value == 0.0 && ignorezeroentries) goto noinsert; 213 if (nonew == 1) goto noinsert; 214 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 215 MatSeqXAIJReallocateAIJ(A,A->rmap.n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 216 N = nrow++ - 1; a->nz++; high++; 217 /* shift up all the later entries in this row */ 218 for (ii=N; ii>=i; ii--) { 219 rp[ii+1] = rp[ii]; 220 ap[ii+1] = ap[ii]; 221 } 222 rp[i] = col; 223 ap[i] = value; 224 noinsert:; 225 low = i + 1; 226 } 227 ailen[row] = nrow; 228 } 229 A->same_nonzero = PETSC_FALSE; 230 PetscFunctionReturn(0); 231 } 232 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "MatGetValues_SeqAIJ" 236 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 237 { 238 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 239 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 240 PetscInt *ai = a->i,*ailen = a->ilen; 241 PetscScalar *ap,*aa = a->a; 242 243 PetscFunctionBegin; 244 for (k=0; k<m; k++) { /* loop over rows */ 245 row = im[k]; 246 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 247 if (row >= A->rmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap.n-1); 248 rp = aj + ai[row]; ap = aa + ai[row]; 249 nrow = ailen[row]; 250 for (l=0; l<n; l++) { /* loop over columns */ 251 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 252 if (in[l] >= A->cmap.n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap.n-1); 253 col = in[l] ; 254 high = nrow; low = 0; /* assume unsorted */ 255 while (high-low > 5) { 256 t = (low+high)/2; 257 if (rp[t] > col) high = t; 258 else low = t; 259 } 260 for (i=low; i<high; i++) { 261 if (rp[i] > col) break; 262 if (rp[i] == col) { 263 *v++ = ap[i]; 264 goto finished; 265 } 266 } 267 *v++ = 0.0; 268 finished:; 269 } 270 } 271 PetscFunctionReturn(0); 272 } 273 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatView_SeqAIJ_Binary" 277 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 278 { 279 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 280 PetscErrorCode ierr; 281 PetscInt i,*col_lens; 282 int fd; 283 284 PetscFunctionBegin; 285 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 286 ierr = PetscMalloc((4+A->rmap.n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 287 col_lens[0] = MAT_FILE_COOKIE; 288 col_lens[1] = A->rmap.n; 289 col_lens[2] = A->cmap.n; 290 col_lens[3] = a->nz; 291 292 /* store lengths of each row and write (including header) to file */ 293 for (i=0; i<A->rmap.n; i++) { 294 col_lens[4+i] = a->i[i+1] - a->i[i]; 295 } 296 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap.n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 297 ierr = PetscFree(col_lens);CHKERRQ(ierr); 298 299 /* store column indices (zero start index) */ 300 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 301 302 /* store nonzero values */ 303 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 308 309 #undef __FUNCT__ 310 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 311 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 312 { 313 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 314 PetscErrorCode ierr; 315 PetscInt i,j,m = A->rmap.n,shift=0; 316 const char *name; 317 PetscViewerFormat format; 318 319 PetscFunctionBegin; 320 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 321 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 322 if (format == PETSC_VIEWER_ASCII_MATLAB) { 323 PetscInt nofinalvalue = 0; 324 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap.n-!shift)) { 325 nofinalvalue = 1; 326 } 327 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap.n);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 330 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 332 333 for (i=0; i<m; i++) { 334 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 335 #if defined(PETSC_USE_COMPLEX) 336 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); 337 #else 338 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 339 #endif 340 } 341 } 342 if (nofinalvalue) { 343 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap.n,0.0);CHKERRQ(ierr); 344 } 345 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 346 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 347 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 348 PetscFunctionReturn(0); 349 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 350 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 351 for (i=0; i<m; i++) { 352 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 353 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 354 #if defined(PETSC_USE_COMPLEX) 355 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 356 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 357 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 358 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 359 } else if (PetscRealPart(a->a[j]) != 0.0) { 360 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 361 } 362 #else 363 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 364 #endif 365 } 366 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 367 } 368 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 369 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 370 PetscInt nzd=0,fshift=1,*sptr; 371 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 372 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 373 for (i=0; i<m; i++) { 374 sptr[i] = nzd+1; 375 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 376 if (a->j[j] >= i) { 377 #if defined(PETSC_USE_COMPLEX) 378 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 379 #else 380 if (a->a[j] != 0.0) nzd++; 381 #endif 382 } 383 } 384 } 385 sptr[m] = nzd+1; 386 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 387 for (i=0; i<m+1; i+=6) { 388 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);} 389 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);} 390 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);} 391 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 392 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 393 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 394 } 395 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 396 ierr = PetscFree(sptr);CHKERRQ(ierr); 397 for (i=0; i<m; i++) { 398 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 399 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 400 } 401 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 402 } 403 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 404 for (i=0; i<m; i++) { 405 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 406 if (a->j[j] >= i) { 407 #if defined(PETSC_USE_COMPLEX) 408 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 409 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 410 } 411 #else 412 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 413 #endif 414 } 415 } 416 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 417 } 418 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 419 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 420 PetscInt cnt = 0,jcnt; 421 PetscScalar value; 422 423 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 424 for (i=0; i<m; i++) { 425 jcnt = 0; 426 for (j=0; j<A->cmap.n; j++) { 427 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 428 value = a->a[cnt++]; 429 jcnt++; 430 } else { 431 value = 0.0; 432 } 433 #if defined(PETSC_USE_COMPLEX) 434 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 435 #else 436 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 437 #endif 438 } 439 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 440 } 441 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 442 } else { 443 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 444 for (i=0; i<m; i++) { 445 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 446 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 447 #if defined(PETSC_USE_COMPLEX) 448 if (PetscImaginaryPart(a->a[j]) > 0.0) { 449 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 450 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 451 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 452 } else { 453 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 454 } 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 457 #endif 458 } 459 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 460 } 461 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 462 } 463 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 464 PetscFunctionReturn(0); 465 } 466 467 #undef __FUNCT__ 468 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 469 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 470 { 471 Mat A = (Mat) Aa; 472 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 473 PetscErrorCode ierr; 474 PetscInt i,j,m = A->rmap.n,color; 475 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 476 PetscViewer viewer; 477 PetscViewerFormat format; 478 479 PetscFunctionBegin; 480 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 481 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 482 483 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 484 /* loop over matrix elements drawing boxes */ 485 486 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 487 /* Blue for negative, Cyan for zero and Red for positive */ 488 color = PETSC_DRAW_BLUE; 489 for (i=0; i<m; i++) { 490 y_l = m - i - 1.0; y_r = y_l + 1.0; 491 for (j=a->i[i]; j<a->i[i+1]; j++) { 492 x_l = a->j[j] ; x_r = x_l + 1.0; 493 #if defined(PETSC_USE_COMPLEX) 494 if (PetscRealPart(a->a[j]) >= 0.) continue; 495 #else 496 if (a->a[j] >= 0.) continue; 497 #endif 498 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 499 } 500 } 501 color = PETSC_DRAW_CYAN; 502 for (i=0; i<m; i++) { 503 y_l = m - i - 1.0; y_r = y_l + 1.0; 504 for (j=a->i[i]; j<a->i[i+1]; j++) { 505 x_l = a->j[j]; x_r = x_l + 1.0; 506 if (a->a[j] != 0.) continue; 507 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 508 } 509 } 510 color = PETSC_DRAW_RED; 511 for (i=0; i<m; i++) { 512 y_l = m - i - 1.0; y_r = y_l + 1.0; 513 for (j=a->i[i]; j<a->i[i+1]; j++) { 514 x_l = a->j[j]; x_r = x_l + 1.0; 515 #if defined(PETSC_USE_COMPLEX) 516 if (PetscRealPart(a->a[j]) <= 0.) continue; 517 #else 518 if (a->a[j] <= 0.) continue; 519 #endif 520 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 521 } 522 } 523 } else { 524 /* use contour shading to indicate magnitude of values */ 525 /* first determine max of all nonzero values */ 526 PetscInt nz = a->nz,count; 527 PetscDraw popup; 528 PetscReal scale; 529 530 for (i=0; i<nz; i++) { 531 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 532 } 533 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 534 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 535 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 536 count = 0; 537 for (i=0; i<m; i++) { 538 y_l = m - i - 1.0; y_r = y_l + 1.0; 539 for (j=a->i[i]; j<a->i[i+1]; j++) { 540 x_l = a->j[j]; x_r = x_l + 1.0; 541 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 542 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 543 count++; 544 } 545 } 546 } 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "MatView_SeqAIJ_Draw" 552 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 553 { 554 PetscErrorCode ierr; 555 PetscDraw draw; 556 PetscReal xr,yr,xl,yl,h,w; 557 PetscTruth isnull; 558 559 PetscFunctionBegin; 560 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 561 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 562 if (isnull) PetscFunctionReturn(0); 563 564 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 565 xr = A->cmap.n; yr = A->rmap.n; h = yr/10.0; w = xr/10.0; 566 xr += w; yr += h; xl = -w; yl = -h; 567 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 568 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 569 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 570 PetscFunctionReturn(0); 571 } 572 573 #undef __FUNCT__ 574 #define __FUNCT__ "MatView_SeqAIJ" 575 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 576 { 577 PetscErrorCode ierr; 578 PetscTruth iascii,isbinary,isdraw; 579 580 PetscFunctionBegin; 581 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 582 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 583 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 584 if (iascii) { 585 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 586 } else if (isbinary) { 587 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 588 } else if (isdraw) { 589 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 590 } else { 591 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 592 } 593 ierr = MatView_Inode(A,viewer);CHKERRQ(ierr); 594 PetscFunctionReturn(0); 595 } 596 597 #undef __FUNCT__ 598 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 599 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 600 { 601 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 602 PetscErrorCode ierr; 603 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 604 PetscInt m = A->rmap.n,*ip,N,*ailen = a->ilen,rmax = 0; 605 PetscScalar *aa = a->a,*ap; 606 PetscReal ratio=0.6; 607 608 PetscFunctionBegin; 609 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 610 611 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 612 for (i=1; i<m; i++) { 613 /* move each row back by the amount of empty slots (fshift) before it*/ 614 fshift += imax[i-1] - ailen[i-1]; 615 rmax = PetscMax(rmax,ailen[i]); 616 if (fshift) { 617 ip = aj + ai[i] ; 618 ap = aa + ai[i] ; 619 N = ailen[i]; 620 for (j=0; j<N; j++) { 621 ip[j-fshift] = ip[j]; 622 ap[j-fshift] = ap[j]; 623 } 624 } 625 ai[i] = ai[i-1] + ailen[i-1]; 626 } 627 if (m) { 628 fshift += imax[m-1] - ailen[m-1]; 629 ai[m] = ai[m-1] + ailen[m-1]; 630 } 631 /* reset ilen and imax for each row */ 632 for (i=0; i<m; i++) { 633 ailen[i] = imax[i] = ai[i+1] - ai[i]; 634 } 635 a->nz = ai[m]; 636 637 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 638 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap.n,fshift,a->nz);CHKERRQ(ierr); 639 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 640 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 641 642 a->reallocs = 0; 643 A->info.nz_unneeded = (double)fshift; 644 a->rmax = rmax; 645 646 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 647 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 648 A->same_nonzero = PETSC_TRUE; 649 650 ierr = MatAssemblyEnd_Inode(A,mode);CHKERRQ(ierr); 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "MatRealPart_SeqAIJ" 656 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 657 { 658 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 659 PetscInt i,nz = a->nz; 660 PetscScalar *aa = a->a; 661 662 PetscFunctionBegin; 663 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 669 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 670 { 671 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 672 PetscInt i,nz = a->nz; 673 PetscScalar *aa = a->a; 674 675 PetscFunctionBegin; 676 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 682 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 683 { 684 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 685 PetscErrorCode ierr; 686 687 PetscFunctionBegin; 688 ierr = PetscMemzero(a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 #undef __FUNCT__ 693 #define __FUNCT__ "MatDestroy_SeqAIJ" 694 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 695 { 696 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 #if defined(PETSC_USE_LOG) 701 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap.n,A->cmap.n,a->nz); 702 #endif 703 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 704 if (a->row) { 705 ierr = ISDestroy(a->row);CHKERRQ(ierr); 706 } 707 if (a->col) { 708 ierr = ISDestroy(a->col);CHKERRQ(ierr); 709 } 710 ierr = PetscFree(a->diag);CHKERRQ(ierr); 711 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 712 ierr = PetscFree(a->idiag);CHKERRQ(ierr); 713 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 714 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 715 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 716 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 717 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 718 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 719 if (a->compressedrow.use){ierr = PetscFree(a->compressedrow.i);} 720 721 ierr = MatDestroy_Inode(A);CHKERRQ(ierr); 722 723 ierr = PetscFree(a);CHKERRQ(ierr); 724 725 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 726 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 727 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 728 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 729 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 730 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 732 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 735 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 736 PetscFunctionReturn(0); 737 } 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "MatCompress_SeqAIJ" 741 PetscErrorCode MatCompress_SeqAIJ(Mat A) 742 { 743 PetscFunctionBegin; 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "MatSetOption_SeqAIJ" 749 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 750 { 751 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 switch (op) { 756 case MAT_ROW_ORIENTED: 757 a->roworiented = flg; 758 break; 759 case MAT_KEEP_ZEROED_ROWS: 760 a->keepzeroedrows = flg; 761 break; 762 case MAT_NEW_NONZERO_LOCATIONS: 763 a->nonew = (flg ? 0 : 1); 764 break; 765 case MAT_NEW_NONZERO_LOCATION_ERR: 766 a->nonew = (flg ? -1 : 0); 767 break; 768 case MAT_NEW_NONZERO_ALLOCATION_ERR: 769 a->nonew = (flg ? -2 : 0); 770 break; 771 case MAT_IGNORE_ZERO_ENTRIES: 772 a->ignorezeroentries = flg; 773 break; 774 case MAT_USE_COMPRESSEDROW: 775 a->compressedrow.use = flg; 776 break; 777 case MAT_NEW_DIAGONALS: 778 case MAT_IGNORE_OFF_PROC_ENTRIES: 779 case MAT_USE_HASH_TABLE: 780 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 781 break; 782 default: 783 break; 784 } 785 ierr = MatSetOption_Inode(A,op,flg);CHKERRQ(ierr); 786 PetscFunctionReturn(0); 787 } 788 789 #undef __FUNCT__ 790 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 791 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 792 { 793 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 794 PetscErrorCode ierr; 795 PetscInt i,j,n; 796 PetscScalar *x,zero = 0.0; 797 798 PetscFunctionBegin; 799 ierr = VecSet(v,zero);CHKERRQ(ierr); 800 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 801 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 802 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 803 for (i=0; i<A->rmap.n; i++) { 804 for (j=a->i[i]; j<a->i[i+1]; j++) { 805 if (a->j[j] == i) { 806 x[i] = a->a[j]; 807 break; 808 } 809 } 810 } 811 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 812 PetscFunctionReturn(0); 813 } 814 815 #undef __FUNCT__ 816 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 817 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 818 { 819 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 820 PetscScalar *x,*y; 821 PetscErrorCode ierr; 822 PetscInt m = A->rmap.n; 823 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 824 PetscScalar *v,alpha; 825 PetscInt n,i,*idx,*ii,*ridx=PETSC_NULL; 826 Mat_CompressedRow cprow = a->compressedrow; 827 PetscTruth usecprow = cprow.use; 828 #endif 829 830 PetscFunctionBegin; 831 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 832 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 833 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 834 835 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 836 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 837 #else 838 if (usecprow){ 839 m = cprow.nrows; 840 ii = cprow.i; 841 ridx = cprow.rindex; 842 } else { 843 ii = a->i; 844 } 845 for (i=0; i<m; i++) { 846 idx = a->j + ii[i] ; 847 v = a->a + ii[i] ; 848 n = ii[i+1] - ii[i]; 849 if (usecprow){ 850 alpha = x[ridx[i]]; 851 } else { 852 alpha = x[i]; 853 } 854 while (n-->0) {y[*idx++] += alpha * *v++;} 855 } 856 #endif 857 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 858 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 859 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 860 PetscFunctionReturn(0); 861 } 862 863 #undef __FUNCT__ 864 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 865 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 866 { 867 PetscScalar zero = 0.0; 868 PetscErrorCode ierr; 869 870 PetscFunctionBegin; 871 ierr = VecSet(yy,zero);CHKERRQ(ierr); 872 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 877 #undef __FUNCT__ 878 #define __FUNCT__ "MatMult_SeqAIJ" 879 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 880 { 881 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 882 PetscScalar *x,*y,*aa; 883 PetscErrorCode ierr; 884 PetscInt m=A->rmap.n,*aj,*ii; 885 PetscInt n,i,j,nonzerorow=0,*ridx=PETSC_NULL; 886 PetscScalar sum; 887 PetscTruth usecprow=a->compressedrow.use; 888 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 889 PetscInt jrow; 890 #endif 891 892 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 893 #pragma disjoint(*x,*y,*aa) 894 #endif 895 896 PetscFunctionBegin; 897 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 898 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 899 aj = a->j; 900 aa = a->a; 901 ii = a->i; 902 if (usecprow){ /* use compressed row format */ 903 m = a->compressedrow.nrows; 904 ii = a->compressedrow.i; 905 ridx = a->compressedrow.rindex; 906 for (i=0; i<m; i++){ 907 n = ii[i+1] - ii[i]; 908 aj = a->j + ii[i]; 909 aa = a->a + ii[i]; 910 sum = 0.0; 911 nonzerorow += (n>0); 912 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 913 y[*ridx++] = sum; 914 } 915 } else { /* do not use compressed row format */ 916 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 917 fortranmultaij_(&m,x,ii,aj,aa,y); 918 #else 919 for (i=0; i<m; i++) { 920 jrow = ii[i]; 921 n = ii[i+1] - jrow; 922 sum = 0.0; 923 nonzerorow += (n>0); 924 for (j=0; j<n; j++) { 925 sum += aa[jrow]*x[aj[jrow]]; jrow++; 926 } 927 y[i] = sum; 928 } 929 #endif 930 } 931 ierr = PetscLogFlops(2*a->nz - nonzerorow);CHKERRQ(ierr); 932 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 933 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 934 PetscFunctionReturn(0); 935 } 936 937 #undef __FUNCT__ 938 #define __FUNCT__ "MatMultAdd_SeqAIJ" 939 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 940 { 941 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 942 PetscScalar *x,*y,*z,*aa; 943 PetscErrorCode ierr; 944 PetscInt m = A->rmap.n,*aj,*ii; 945 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 946 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 947 PetscScalar sum; 948 PetscTruth usecprow=a->compressedrow.use; 949 #endif 950 951 PetscFunctionBegin; 952 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 953 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 954 if (zz != yy) { 955 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 956 } else { 957 z = y; 958 } 959 960 aj = a->j; 961 aa = a->a; 962 ii = a->i; 963 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 964 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 965 #else 966 if (usecprow){ /* use compressed row format */ 967 if (zz != yy){ 968 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 969 } 970 m = a->compressedrow.nrows; 971 ii = a->compressedrow.i; 972 ridx = a->compressedrow.rindex; 973 for (i=0; i<m; i++){ 974 n = ii[i+1] - ii[i]; 975 aj = a->j + ii[i]; 976 aa = a->a + ii[i]; 977 sum = y[*ridx]; 978 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 979 z[*ridx++] = sum; 980 } 981 } else { /* do not use compressed row format */ 982 for (i=0; i<m; i++) { 983 jrow = ii[i]; 984 n = ii[i+1] - jrow; 985 sum = y[i]; 986 for (j=0; j<n; j++) { 987 sum += aa[jrow]*x[aj[jrow]]; jrow++; 988 } 989 z[i] = sum; 990 } 991 } 992 #endif 993 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 994 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 995 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 996 if (zz != yy) { 997 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 998 } 999 PetscFunctionReturn(0); 1000 } 1001 1002 /* 1003 Adds diagonal pointers to sparse matrix structure. 1004 */ 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1007 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1008 { 1009 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1010 PetscErrorCode ierr; 1011 PetscInt i,j,m = A->rmap.n; 1012 1013 PetscFunctionBegin; 1014 if (!a->diag) { 1015 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1016 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1017 } 1018 for (i=0; i<A->rmap.n; i++) { 1019 a->diag[i] = a->i[i+1]; 1020 for (j=a->i[i]; j<a->i[i+1]; j++) { 1021 if (a->j[j] == i) { 1022 a->diag[i] = j; 1023 break; 1024 } 1025 } 1026 } 1027 PetscFunctionReturn(0); 1028 } 1029 1030 /* 1031 Checks for missing diagonals 1032 */ 1033 #undef __FUNCT__ 1034 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1035 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1036 { 1037 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1038 PetscInt *diag,*jj = a->j,i; 1039 1040 PetscFunctionBegin; 1041 *missing = PETSC_FALSE; 1042 if (A->rmap.n > 0 && !jj) { 1043 *missing = PETSC_TRUE; 1044 if (d) *d = 0; 1045 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1046 } else { 1047 diag = a->diag; 1048 for (i=0; i<A->rmap.n; i++) { 1049 if (jj[diag[i]] != i) { 1050 *missing = PETSC_TRUE; 1051 if (d) *d = i; 1052 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1053 } 1054 } 1055 } 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "MatRelax_SeqAIJ" 1061 PetscErrorCode MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1062 { 1063 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1064 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1065 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1066 PetscErrorCode ierr; 1067 PetscInt n = A->cmap.n,m = A->rmap.n,i; 1068 const PetscInt *idx,*diag; 1069 1070 PetscFunctionBegin; 1071 its = its*lits; 1072 1073 diag = a->diag; 1074 if (!a->idiag) { 1075 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1076 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1077 a->ssor = a->idiag + m; 1078 mdiag = a->ssor + m; 1079 v = a->a; 1080 1081 /* this is wrong when fshift omega changes each iteration */ 1082 if (omega == 1.0 && !fshift) { 1083 for (i=0; i<m; i++) { 1084 mdiag[i] = v[diag[i]]; 1085 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1086 a->idiag[i] = 1.0/v[diag[i]]; 1087 } 1088 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1089 } else { 1090 for (i=0; i<m; i++) { 1091 mdiag[i] = v[diag[i]]; 1092 a->idiag[i] = omega/(fshift + v[diag[i]]); 1093 } 1094 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1095 } 1096 } 1097 t = a->ssor; 1098 idiag = a->idiag; 1099 mdiag = a->idiag + 2*m; 1100 1101 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1102 if (xx != bb) { 1103 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1104 } else { 1105 b = x; 1106 } 1107 1108 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1109 xs = x; 1110 if (flag == SOR_APPLY_UPPER) { 1111 /* apply (U + D/omega) to the vector */ 1112 bs = b; 1113 for (i=0; i<m; i++) { 1114 d = fshift + a->a[diag[i]]; 1115 n = a->i[i+1] - diag[i] - 1; 1116 idx = a->j + diag[i] + 1; 1117 v = a->a + diag[i] + 1; 1118 sum = b[i]*d/omega; 1119 SPARSEDENSEDOT(sum,bs,v,idx,n); 1120 x[i] = sum; 1121 } 1122 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1123 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1124 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 1129 /* Let A = L + U + D; where L is lower trianglar, 1130 U is upper triangular, E is diagonal; This routine applies 1131 1132 (L + E)^{-1} A (U + E)^{-1} 1133 1134 to a vector efficiently using Eisenstat's trick. This is for 1135 the case of SSOR preconditioner, so E is D/omega where omega 1136 is the relaxation factor. 1137 */ 1138 1139 if (flag == SOR_APPLY_LOWER) { 1140 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1141 } else if (flag & SOR_EISENSTAT) { 1142 /* Let A = L + U + D; where L is lower trianglar, 1143 U is upper triangular, E is diagonal; This routine applies 1144 1145 (L + E)^{-1} A (U + E)^{-1} 1146 1147 to a vector efficiently using Eisenstat's trick. This is for 1148 the case of SSOR preconditioner, so E is D/omega where omega 1149 is the relaxation factor. 1150 */ 1151 scale = (2.0/omega) - 1.0; 1152 1153 /* x = (E + U)^{-1} b */ 1154 for (i=m-1; i>=0; i--) { 1155 n = a->i[i+1] - diag[i] - 1; 1156 idx = a->j + diag[i] + 1; 1157 v = a->a + diag[i] + 1; 1158 sum = b[i]; 1159 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1160 x[i] = sum*idiag[i]; 1161 } 1162 1163 /* t = b - (2*E - D)x */ 1164 v = a->a; 1165 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1166 1167 /* t = (E + L)^{-1}t */ 1168 ts = t; 1169 diag = a->diag; 1170 for (i=0; i<m; i++) { 1171 n = diag[i] - a->i[i]; 1172 idx = a->j + a->i[i]; 1173 v = a->a + a->i[i]; 1174 sum = t[i]; 1175 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1176 t[i] = sum*idiag[i]; 1177 /* x = x + t */ 1178 x[i] += t[i]; 1179 } 1180 1181 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1182 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1183 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1184 PetscFunctionReturn(0); 1185 } 1186 if (flag & SOR_ZERO_INITIAL_GUESS) { 1187 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1188 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1189 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1190 #else 1191 for (i=0; i<m; i++) { 1192 n = diag[i] - a->i[i]; 1193 idx = a->j + a->i[i]; 1194 v = a->a + a->i[i]; 1195 sum = b[i]; 1196 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1197 x[i] = sum*idiag[i]; 1198 } 1199 #endif 1200 xb = x; 1201 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1202 } else xb = b; 1203 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1204 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1205 for (i=0; i<m; i++) { 1206 x[i] *= mdiag[i]; 1207 } 1208 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1209 } 1210 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1211 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1212 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1213 #else 1214 for (i=m-1; i>=0; i--) { 1215 n = a->i[i+1] - diag[i] - 1; 1216 idx = a->j + diag[i] + 1; 1217 v = a->a + diag[i] + 1; 1218 sum = xb[i]; 1219 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1220 x[i] = sum*idiag[i]; 1221 } 1222 #endif 1223 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1224 } 1225 its--; 1226 } 1227 while (its--) { 1228 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1229 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1230 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1231 #else 1232 for (i=0; i<m; i++) { 1233 n = a->i[i+1] - a->i[i]; 1234 idx = a->j + a->i[i]; 1235 v = a->a + a->i[i]; 1236 sum = b[i]; 1237 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1238 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1239 } 1240 #endif 1241 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1242 } 1243 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1244 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1245 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1246 #else 1247 for (i=m-1; i>=0; i--) { 1248 n = a->i[i+1] - a->i[i]; 1249 idx = a->j + a->i[i]; 1250 v = a->a + a->i[i]; 1251 sum = b[i]; 1252 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1253 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1254 } 1255 #endif 1256 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1257 } 1258 } 1259 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1260 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1261 PetscFunctionReturn(0); 1262 } 1263 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1267 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1268 { 1269 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1270 1271 PetscFunctionBegin; 1272 info->rows_global = (double)A->rmap.n; 1273 info->columns_global = (double)A->cmap.n; 1274 info->rows_local = (double)A->rmap.n; 1275 info->columns_local = (double)A->cmap.n; 1276 info->block_size = 1.0; 1277 info->nz_allocated = (double)a->maxnz; 1278 info->nz_used = (double)a->nz; 1279 info->nz_unneeded = (double)(a->maxnz - a->nz); 1280 info->assemblies = (double)A->num_ass; 1281 info->mallocs = (double)a->reallocs; 1282 info->memory = A->mem; 1283 if (A->factor) { 1284 info->fill_ratio_given = A->info.fill_ratio_given; 1285 info->fill_ratio_needed = A->info.fill_ratio_needed; 1286 info->factor_mallocs = A->info.factor_mallocs; 1287 } else { 1288 info->fill_ratio_given = 0; 1289 info->fill_ratio_needed = 0; 1290 info->factor_mallocs = 0; 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1297 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1298 { 1299 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1300 PetscInt i,m = A->rmap.n - 1,d; 1301 PetscErrorCode ierr; 1302 PetscTruth missing; 1303 1304 PetscFunctionBegin; 1305 if (a->keepzeroedrows) { 1306 for (i=0; i<N; i++) { 1307 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1308 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1309 } 1310 if (diag != 0.0) { 1311 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1312 if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1313 for (i=0; i<N; i++) { 1314 a->a[a->diag[rows[i]]] = diag; 1315 } 1316 } 1317 A->same_nonzero = PETSC_TRUE; 1318 } else { 1319 if (diag != 0.0) { 1320 for (i=0; i<N; i++) { 1321 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1322 if (a->ilen[rows[i]] > 0) { 1323 a->ilen[rows[i]] = 1; 1324 a->a[a->i[rows[i]]] = diag; 1325 a->j[a->i[rows[i]]] = rows[i]; 1326 } else { /* in case row was completely empty */ 1327 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1328 } 1329 } 1330 } else { 1331 for (i=0; i<N; i++) { 1332 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1333 a->ilen[rows[i]] = 0; 1334 } 1335 } 1336 A->same_nonzero = PETSC_FALSE; 1337 } 1338 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "MatGetRow_SeqAIJ" 1344 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1345 { 1346 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1347 PetscInt *itmp; 1348 1349 PetscFunctionBegin; 1350 if (row < 0 || row >= A->rmap.n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1351 1352 *nz = a->i[row+1] - a->i[row]; 1353 if (v) *v = a->a + a->i[row]; 1354 if (idx) { 1355 itmp = a->j + a->i[row]; 1356 if (*nz) { 1357 *idx = itmp; 1358 } 1359 else *idx = 0; 1360 } 1361 PetscFunctionReturn(0); 1362 } 1363 1364 /* remove this function? */ 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1367 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1368 { 1369 PetscFunctionBegin; 1370 PetscFunctionReturn(0); 1371 } 1372 1373 #undef __FUNCT__ 1374 #define __FUNCT__ "MatNorm_SeqAIJ" 1375 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1376 { 1377 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1378 PetscScalar *v = a->a; 1379 PetscReal sum = 0.0; 1380 PetscErrorCode ierr; 1381 PetscInt i,j; 1382 1383 PetscFunctionBegin; 1384 if (type == NORM_FROBENIUS) { 1385 for (i=0; i<a->nz; i++) { 1386 #if defined(PETSC_USE_COMPLEX) 1387 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1388 #else 1389 sum += (*v)*(*v); v++; 1390 #endif 1391 } 1392 *nrm = sqrt(sum); 1393 } else if (type == NORM_1) { 1394 PetscReal *tmp; 1395 PetscInt *jj = a->j; 1396 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1397 ierr = PetscMemzero(tmp,A->cmap.n*sizeof(PetscReal));CHKERRQ(ierr); 1398 *nrm = 0.0; 1399 for (j=0; j<a->nz; j++) { 1400 tmp[*jj++] += PetscAbsScalar(*v); v++; 1401 } 1402 for (j=0; j<A->cmap.n; j++) { 1403 if (tmp[j] > *nrm) *nrm = tmp[j]; 1404 } 1405 ierr = PetscFree(tmp);CHKERRQ(ierr); 1406 } else if (type == NORM_INFINITY) { 1407 *nrm = 0.0; 1408 for (j=0; j<A->rmap.n; j++) { 1409 v = a->a + a->i[j]; 1410 sum = 0.0; 1411 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1412 sum += PetscAbsScalar(*v); v++; 1413 } 1414 if (sum > *nrm) *nrm = sum; 1415 } 1416 } else { 1417 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1418 } 1419 PetscFunctionReturn(0); 1420 } 1421 1422 #undef __FUNCT__ 1423 #define __FUNCT__ "MatTranspose_SeqAIJ" 1424 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1425 { 1426 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1427 Mat C; 1428 PetscErrorCode ierr; 1429 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap.n,len,*col; 1430 PetscScalar *array = a->a; 1431 1432 PetscFunctionBegin; 1433 if (!B && m != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1434 ierr = PetscMalloc((1+A->cmap.n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1435 ierr = PetscMemzero(col,(1+A->cmap.n)*sizeof(PetscInt));CHKERRQ(ierr); 1436 1437 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1438 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1439 ierr = MatSetSizes(C,A->cmap.n,m,A->cmap.n,m);CHKERRQ(ierr); 1440 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1441 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1442 ierr = PetscFree(col);CHKERRQ(ierr); 1443 for (i=0; i<m; i++) { 1444 len = ai[i+1]-ai[i]; 1445 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1446 array += len; 1447 aj += len; 1448 } 1449 1450 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1451 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1452 1453 if (B) { 1454 *B = C; 1455 } else { 1456 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1457 } 1458 PetscFunctionReturn(0); 1459 } 1460 1461 EXTERN_C_BEGIN 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1464 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1465 { 1466 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1467 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1468 PetscErrorCode ierr; 1469 PetscInt ma,na,mb,nb, i; 1470 1471 PetscFunctionBegin; 1472 bij = (Mat_SeqAIJ *) B->data; 1473 1474 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1475 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1476 if (ma!=nb || na!=mb){ 1477 *f = PETSC_FALSE; 1478 PetscFunctionReturn(0); 1479 } 1480 aii = aij->i; bii = bij->i; 1481 adx = aij->j; bdx = bij->j; 1482 va = aij->a; vb = bij->a; 1483 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1484 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1485 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1486 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1487 1488 *f = PETSC_TRUE; 1489 for (i=0; i<ma; i++) { 1490 while (aptr[i]<aii[i+1]) { 1491 PetscInt idc,idr; 1492 PetscScalar vc,vr; 1493 /* column/row index/value */ 1494 idc = adx[aptr[i]]; 1495 idr = bdx[bptr[idc]]; 1496 vc = va[aptr[i]]; 1497 vr = vb[bptr[idc]]; 1498 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1499 *f = PETSC_FALSE; 1500 goto done; 1501 } else { 1502 aptr[i]++; 1503 if (B || i!=idc) bptr[idc]++; 1504 } 1505 } 1506 } 1507 done: 1508 ierr = PetscFree(aptr);CHKERRQ(ierr); 1509 if (B) { 1510 ierr = PetscFree(bptr);CHKERRQ(ierr); 1511 } 1512 PetscFunctionReturn(0); 1513 } 1514 EXTERN_C_END 1515 1516 EXTERN_C_BEGIN 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1519 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1520 { 1521 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1522 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1523 PetscErrorCode ierr; 1524 PetscInt ma,na,mb,nb, i; 1525 1526 PetscFunctionBegin; 1527 bij = (Mat_SeqAIJ *) B->data; 1528 1529 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1530 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1531 if (ma!=nb || na!=mb){ 1532 *f = PETSC_FALSE; 1533 PetscFunctionReturn(0); 1534 } 1535 aii = aij->i; bii = bij->i; 1536 adx = aij->j; bdx = bij->j; 1537 va = aij->a; vb = bij->a; 1538 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1539 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1540 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1541 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1542 1543 *f = PETSC_TRUE; 1544 for (i=0; i<ma; i++) { 1545 while (aptr[i]<aii[i+1]) { 1546 PetscInt idc,idr; 1547 PetscScalar vc,vr; 1548 /* column/row index/value */ 1549 idc = adx[aptr[i]]; 1550 idr = bdx[bptr[idc]]; 1551 vc = va[aptr[i]]; 1552 vr = vb[bptr[idc]]; 1553 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1554 *f = PETSC_FALSE; 1555 goto done; 1556 } else { 1557 aptr[i]++; 1558 if (B || i!=idc) bptr[idc]++; 1559 } 1560 } 1561 } 1562 done: 1563 ierr = PetscFree(aptr);CHKERRQ(ierr); 1564 if (B) { 1565 ierr = PetscFree(bptr);CHKERRQ(ierr); 1566 } 1567 PetscFunctionReturn(0); 1568 } 1569 EXTERN_C_END 1570 1571 #undef __FUNCT__ 1572 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1573 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1574 { 1575 PetscErrorCode ierr; 1576 PetscFunctionBegin; 1577 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1583 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1584 { 1585 PetscErrorCode ierr; 1586 PetscFunctionBegin; 1587 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 1591 #undef __FUNCT__ 1592 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1593 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1594 { 1595 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1596 PetscScalar *l,*r,x,*v; 1597 PetscErrorCode ierr; 1598 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,M,nz = a->nz,*jj; 1599 1600 PetscFunctionBegin; 1601 if (ll) { 1602 /* The local size is used so that VecMPI can be passed to this routine 1603 by MatDiagonalScale_MPIAIJ */ 1604 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1605 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1606 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1607 v = a->a; 1608 for (i=0; i<m; i++) { 1609 x = l[i]; 1610 M = a->i[i+1] - a->i[i]; 1611 for (j=0; j<M; j++) { (*v++) *= x;} 1612 } 1613 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1614 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1615 } 1616 if (rr) { 1617 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1618 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1619 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1620 v = a->a; jj = a->j; 1621 for (i=0; i<nz; i++) { 1622 (*v++) *= r[*jj++]; 1623 } 1624 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1625 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1626 } 1627 PetscFunctionReturn(0); 1628 } 1629 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1632 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1633 { 1634 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1635 PetscErrorCode ierr; 1636 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap.n,*lens; 1637 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1638 PetscInt *irow,*icol,nrows,ncols; 1639 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1640 PetscScalar *a_new,*mat_a; 1641 Mat C; 1642 PetscTruth stride; 1643 1644 PetscFunctionBegin; 1645 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1646 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1647 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1648 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1649 1650 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1651 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1652 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1653 1654 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1655 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1656 if (stride && step == 1) { 1657 /* special case of contiguous rows */ 1658 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1659 starts = lens + nrows; 1660 /* loop over new rows determining lens and starting points */ 1661 for (i=0; i<nrows; i++) { 1662 kstart = ai[irow[i]]; 1663 kend = kstart + ailen[irow[i]]; 1664 for (k=kstart; k<kend; k++) { 1665 if (aj[k] >= first) { 1666 starts[i] = k; 1667 break; 1668 } 1669 } 1670 sum = 0; 1671 while (k < kend) { 1672 if (aj[k++] >= first+ncols) break; 1673 sum++; 1674 } 1675 lens[i] = sum; 1676 } 1677 /* create submatrix */ 1678 if (scall == MAT_REUSE_MATRIX) { 1679 PetscInt n_cols,n_rows; 1680 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1681 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1682 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1683 C = *B; 1684 } else { 1685 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1686 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1687 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1688 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1689 } 1690 c = (Mat_SeqAIJ*)C->data; 1691 1692 /* loop over rows inserting into submatrix */ 1693 a_new = c->a; 1694 j_new = c->j; 1695 i_new = c->i; 1696 1697 for (i=0; i<nrows; i++) { 1698 ii = starts[i]; 1699 lensi = lens[i]; 1700 for (k=0; k<lensi; k++) { 1701 *j_new++ = aj[ii+k] - first; 1702 } 1703 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1704 a_new += lensi; 1705 i_new[i+1] = i_new[i] + lensi; 1706 c->ilen[i] = lensi; 1707 } 1708 ierr = PetscFree(lens);CHKERRQ(ierr); 1709 } else { 1710 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1711 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1712 1713 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1714 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1715 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1716 /* determine lens of each row */ 1717 for (i=0; i<nrows; i++) { 1718 kstart = ai[irow[i]]; 1719 kend = kstart + a->ilen[irow[i]]; 1720 lens[i] = 0; 1721 for (k=kstart; k<kend; k++) { 1722 if (smap[aj[k]]) { 1723 lens[i]++; 1724 } 1725 } 1726 } 1727 /* Create and fill new matrix */ 1728 if (scall == MAT_REUSE_MATRIX) { 1729 PetscTruth equal; 1730 1731 c = (Mat_SeqAIJ *)((*B)->data); 1732 if ((*B)->rmap.n != nrows || (*B)->cmap.n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1733 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap.n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1734 if (!equal) { 1735 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1736 } 1737 ierr = PetscMemzero(c->ilen,(*B)->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 1738 C = *B; 1739 } else { 1740 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1741 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1742 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1743 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1744 } 1745 c = (Mat_SeqAIJ *)(C->data); 1746 for (i=0; i<nrows; i++) { 1747 row = irow[i]; 1748 kstart = ai[row]; 1749 kend = kstart + a->ilen[row]; 1750 mat_i = c->i[i]; 1751 mat_j = c->j + mat_i; 1752 mat_a = c->a + mat_i; 1753 mat_ilen = c->ilen + i; 1754 for (k=kstart; k<kend; k++) { 1755 if ((tcol=smap[a->j[k]])) { 1756 *mat_j++ = tcol - 1; 1757 *mat_a++ = a->a[k]; 1758 (*mat_ilen)++; 1759 1760 } 1761 } 1762 } 1763 /* Free work space */ 1764 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1765 ierr = PetscFree(smap);CHKERRQ(ierr); 1766 ierr = PetscFree(lens);CHKERRQ(ierr); 1767 } 1768 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1769 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1770 1771 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1772 *B = C; 1773 PetscFunctionReturn(0); 1774 } 1775 1776 /* 1777 */ 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1780 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1781 { 1782 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1783 PetscErrorCode ierr; 1784 Mat outA; 1785 PetscTruth row_identity,col_identity; 1786 1787 PetscFunctionBegin; 1788 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1789 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1790 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1791 1792 outA = inA; 1793 inA->factor = FACTOR_LU; 1794 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1795 if (a->row) { ierr = ISDestroy(a->row); CHKERRQ(ierr);} 1796 a->row = row; 1797 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1798 if (a->col) { ierr = ISDestroy(a->col); CHKERRQ(ierr);} 1799 a->col = col; 1800 1801 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1802 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1803 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1804 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1805 1806 if (!a->solve_work) { /* this matrix may have been factored before */ 1807 ierr = PetscMalloc((inA->rmap.n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1808 ierr = PetscLogObjectMemory(inA, (inA->rmap.n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1809 } 1810 1811 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1812 if (row_identity && col_identity) { 1813 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1814 } else { 1815 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(inA,info,&outA);CHKERRQ(ierr); 1816 } 1817 PetscFunctionReturn(0); 1818 } 1819 1820 #include "petscblaslapack.h" 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "MatScale_SeqAIJ" 1823 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1824 { 1825 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1826 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1827 PetscScalar oalpha = alpha; 1828 PetscErrorCode ierr; 1829 1830 1831 PetscFunctionBegin; 1832 BLASscal_(&bnz,&oalpha,a->a,&one); 1833 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1834 PetscFunctionReturn(0); 1835 } 1836 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1839 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1840 { 1841 PetscErrorCode ierr; 1842 PetscInt i; 1843 1844 PetscFunctionBegin; 1845 if (scall == MAT_INITIAL_MATRIX) { 1846 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1847 } 1848 1849 for (i=0; i<n; i++) { 1850 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1851 } 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1857 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1858 { 1859 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1860 PetscErrorCode ierr; 1861 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1862 PetscInt start,end,*ai,*aj; 1863 PetscBT table; 1864 1865 PetscFunctionBegin; 1866 m = A->rmap.n; 1867 ai = a->i; 1868 aj = a->j; 1869 1870 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1871 1872 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1873 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1874 1875 for (i=0; i<is_max; i++) { 1876 /* Initialize the two local arrays */ 1877 isz = 0; 1878 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1879 1880 /* Extract the indices, assume there can be duplicate entries */ 1881 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1882 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1883 1884 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1885 for (j=0; j<n ; ++j){ 1886 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1887 } 1888 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1889 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1890 1891 k = 0; 1892 for (j=0; j<ov; j++){ /* for each overlap */ 1893 n = isz; 1894 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1895 row = nidx[k]; 1896 start = ai[row]; 1897 end = ai[row+1]; 1898 for (l = start; l<end ; l++){ 1899 val = aj[l] ; 1900 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1901 } 1902 } 1903 } 1904 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1905 } 1906 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1907 ierr = PetscFree(nidx);CHKERRQ(ierr); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /* -------------------------------------------------------------- */ 1912 #undef __FUNCT__ 1913 #define __FUNCT__ "MatPermute_SeqAIJ" 1914 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1915 { 1916 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1917 PetscErrorCode ierr; 1918 PetscInt i,nz,m = A->rmap.n,n = A->cmap.n,*col; 1919 PetscInt *row,*cnew,j,*lens; 1920 IS icolp,irowp; 1921 PetscInt *cwork; 1922 PetscScalar *vwork; 1923 1924 PetscFunctionBegin; 1925 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1926 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1927 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1928 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1929 1930 /* determine lengths of permuted rows */ 1931 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1932 for (i=0; i<m; i++) { 1933 lens[row[i]] = a->i[i+1] - a->i[i]; 1934 } 1935 ierr = MatCreate(A->comm,B);CHKERRQ(ierr); 1936 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1937 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1938 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1939 ierr = PetscFree(lens);CHKERRQ(ierr); 1940 1941 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1942 for (i=0; i<m; i++) { 1943 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1944 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1945 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1946 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1947 } 1948 ierr = PetscFree(cnew);CHKERRQ(ierr); 1949 (*B)->assembled = PETSC_FALSE; 1950 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1951 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1952 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1953 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1954 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1955 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 #undef __FUNCT__ 1960 #define __FUNCT__ "MatCopy_SeqAIJ" 1961 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1962 { 1963 PetscErrorCode ierr; 1964 1965 PetscFunctionBegin; 1966 /* If the two matrices have the same copy implementation, use fast copy. */ 1967 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1968 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1969 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1970 1971 if (a->i[A->rmap.n] != b->i[B->rmap.n]) { 1972 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1973 } 1974 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 1975 } else { 1976 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1977 } 1978 PetscFunctionReturn(0); 1979 } 1980 1981 #undef __FUNCT__ 1982 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1983 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1984 { 1985 PetscErrorCode ierr; 1986 1987 PetscFunctionBegin; 1988 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1989 PetscFunctionReturn(0); 1990 } 1991 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "MatGetArray_SeqAIJ" 1994 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1995 { 1996 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1997 PetscFunctionBegin; 1998 *array = a->a; 1999 PetscFunctionReturn(0); 2000 } 2001 2002 #undef __FUNCT__ 2003 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2004 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2005 { 2006 PetscFunctionBegin; 2007 PetscFunctionReturn(0); 2008 } 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2012 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2013 { 2014 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2015 PetscErrorCode ierr; 2016 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2017 PetscScalar dx,*y,*xx,*w3_array; 2018 PetscScalar *vscale_array; 2019 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2020 Vec w1,w2,w3; 2021 void *fctx = coloring->fctx; 2022 PetscTruth flg; 2023 2024 PetscFunctionBegin; 2025 if (!coloring->w1) { 2026 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2027 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2028 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2029 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2030 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2031 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2032 } 2033 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2034 2035 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2036 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 2037 if (flg) { 2038 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2039 } else { 2040 PetscTruth assembled; 2041 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2042 if (assembled) { 2043 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2044 } 2045 } 2046 2047 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2048 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2049 2050 /* 2051 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2052 coloring->F for the coarser grids from the finest 2053 */ 2054 if (coloring->F) { 2055 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2056 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2057 if (m1 != m2) { 2058 coloring->F = 0; 2059 } 2060 } 2061 2062 if (coloring->F) { 2063 w1 = coloring->F; 2064 coloring->F = 0; 2065 } else { 2066 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2067 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2068 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2069 } 2070 2071 /* 2072 Compute all the scale factors and share with other processors 2073 */ 2074 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2075 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2076 for (k=0; k<coloring->ncolors; k++) { 2077 /* 2078 Loop over each column associated with color adding the 2079 perturbation to the vector w3. 2080 */ 2081 for (l=0; l<coloring->ncolumns[k]; l++) { 2082 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2083 dx = xx[col]; 2084 if (dx == 0.0) dx = 1.0; 2085 #if !defined(PETSC_USE_COMPLEX) 2086 if (dx < umin && dx >= 0.0) dx = umin; 2087 else if (dx < 0.0 && dx > -umin) dx = -umin; 2088 #else 2089 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2090 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2091 #endif 2092 dx *= epsilon; 2093 vscale_array[col] = 1.0/dx; 2094 } 2095 } 2096 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2097 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2098 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2099 2100 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2101 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2102 2103 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2104 else vscaleforrow = coloring->columnsforrow; 2105 2106 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2107 /* 2108 Loop over each color 2109 */ 2110 for (k=0; k<coloring->ncolors; k++) { 2111 coloring->currentcolor = k; 2112 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2113 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2114 /* 2115 Loop over each column associated with color adding the 2116 perturbation to the vector w3. 2117 */ 2118 for (l=0; l<coloring->ncolumns[k]; l++) { 2119 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2120 dx = xx[col]; 2121 if (dx == 0.0) dx = 1.0; 2122 #if !defined(PETSC_USE_COMPLEX) 2123 if (dx < umin && dx >= 0.0) dx = umin; 2124 else if (dx < 0.0 && dx > -umin) dx = -umin; 2125 #else 2126 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2127 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2128 #endif 2129 dx *= epsilon; 2130 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2131 w3_array[col] += dx; 2132 } 2133 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2134 2135 /* 2136 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2137 */ 2138 2139 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2140 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2141 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2142 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2143 2144 /* 2145 Loop over rows of vector, putting results into Jacobian matrix 2146 */ 2147 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2148 for (l=0; l<coloring->nrows[k]; l++) { 2149 row = coloring->rows[k][l]; 2150 col = coloring->columnsforrow[k][l]; 2151 y[row] *= vscale_array[vscaleforrow[k][l]]; 2152 srow = row + start; 2153 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2154 } 2155 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2156 } 2157 coloring->currentcolor = k; 2158 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2159 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2160 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2161 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2162 PetscFunctionReturn(0); 2163 } 2164 2165 #include "petscblaslapack.h" 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "MatAXPY_SeqAIJ" 2168 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2169 { 2170 PetscErrorCode ierr; 2171 PetscInt i; 2172 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2173 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2174 2175 PetscFunctionBegin; 2176 if (str == SAME_NONZERO_PATTERN) { 2177 PetscScalar alpha = a; 2178 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2179 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2180 if (y->xtoy && y->XtoY != X) { 2181 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2182 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2183 } 2184 if (!y->xtoy) { /* get xtoy */ 2185 ierr = MatAXPYGetxtoy_Private(X->rmap.n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2186 y->XtoY = X; 2187 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2188 } 2189 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2190 ierr = PetscInfo3(0,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2191 } else { 2192 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2193 } 2194 PetscFunctionReturn(0); 2195 } 2196 2197 #undef __FUNCT__ 2198 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2199 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2200 { 2201 PetscFunctionBegin; 2202 PetscFunctionReturn(0); 2203 } 2204 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "MatConjugate_SeqAIJ" 2207 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2208 { 2209 #if defined(PETSC_USE_COMPLEX) 2210 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2211 PetscInt i,nz; 2212 PetscScalar *a; 2213 2214 PetscFunctionBegin; 2215 nz = aij->nz; 2216 a = aij->a; 2217 for (i=0; i<nz; i++) { 2218 a[i] = PetscConj(a[i]); 2219 } 2220 #else 2221 PetscFunctionBegin; 2222 #endif 2223 PetscFunctionReturn(0); 2224 } 2225 2226 #undef __FUNCT__ 2227 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2228 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2229 { 2230 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2231 PetscErrorCode ierr; 2232 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2233 PetscReal atmp; 2234 PetscScalar *x; 2235 MatScalar *aa; 2236 2237 PetscFunctionBegin; 2238 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2239 aa = a->a; 2240 ai = a->i; 2241 aj = a->j; 2242 2243 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2244 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2245 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2246 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2247 for (i=0; i<m; i++) { 2248 ncols = ai[1] - ai[0]; ai++; 2249 x[i] = 0.0; if (idx) idx[i] = 0; 2250 for (j=0; j<ncols; j++){ 2251 atmp = PetscAbsScalar(*aa); 2252 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2253 aa++; aj++; 2254 } 2255 } 2256 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 #undef __FUNCT__ 2261 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2262 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2263 { 2264 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2265 PetscErrorCode ierr; 2266 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2267 PetscScalar *x; 2268 MatScalar *aa; 2269 2270 PetscFunctionBegin; 2271 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2272 aa = a->a; 2273 ai = a->i; 2274 aj = a->j; 2275 2276 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2277 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2278 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2279 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2280 for (i=0; i<m; i++) { 2281 ncols = ai[1] - ai[0]; ai++; 2282 if (ncols == A->cmap.n) { /* row is dense */ 2283 x[i] = *aa; if (idx) idx[i] = 0; 2284 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2285 x[i] = 0.0; 2286 if (idx) { 2287 idx[i] = 0; /* in case ncols is zero */ 2288 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2289 if (aj[j] > j) { 2290 idx[i] = j; 2291 break; 2292 } 2293 } 2294 } 2295 } 2296 for (j=0; j<ncols; j++){ 2297 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2298 aa++; aj++; 2299 } 2300 } 2301 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2302 PetscFunctionReturn(0); 2303 } 2304 2305 #undef __FUNCT__ 2306 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2307 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2308 { 2309 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2310 PetscErrorCode ierr; 2311 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2312 PetscScalar *x; 2313 MatScalar *aa; 2314 2315 PetscFunctionBegin; 2316 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2317 aa = a->a; 2318 ai = a->i; 2319 aj = a->j; 2320 2321 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2322 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2323 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2324 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2325 for (i=0; i<m; i++) { 2326 ncols = ai[1] - ai[0]; ai++; 2327 if (ncols == A->cmap.n) { /* row is dense */ 2328 x[i] = *aa; if (idx) idx[i] = 0; 2329 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2330 x[i] = 0.0; 2331 if (idx) { /* find first implicit 0.0 in the row */ 2332 idx[i] = 0; /* in case ncols is zero */ 2333 for (j=0;j<ncols;j++) { 2334 if (aj[j] > j) { 2335 idx[i] = j; 2336 break; 2337 } 2338 } 2339 } 2340 } 2341 for (j=0; j<ncols; j++){ 2342 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2343 aa++; aj++; 2344 } 2345 } 2346 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2347 PetscFunctionReturn(0); 2348 } 2349 2350 /* -------------------------------------------------------------------*/ 2351 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2352 MatGetRow_SeqAIJ, 2353 MatRestoreRow_SeqAIJ, 2354 MatMult_SeqAIJ, 2355 /* 4*/ MatMultAdd_SeqAIJ, 2356 MatMultTranspose_SeqAIJ, 2357 MatMultTransposeAdd_SeqAIJ, 2358 MatSolve_SeqAIJ, 2359 MatSolveAdd_SeqAIJ, 2360 MatSolveTranspose_SeqAIJ, 2361 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2362 MatLUFactor_SeqAIJ, 2363 0, 2364 MatRelax_SeqAIJ, 2365 MatTranspose_SeqAIJ, 2366 /*15*/ MatGetInfo_SeqAIJ, 2367 MatEqual_SeqAIJ, 2368 MatGetDiagonal_SeqAIJ, 2369 MatDiagonalScale_SeqAIJ, 2370 MatNorm_SeqAIJ, 2371 /*20*/ 0, 2372 MatAssemblyEnd_SeqAIJ, 2373 MatCompress_SeqAIJ, 2374 MatSetOption_SeqAIJ, 2375 MatZeroEntries_SeqAIJ, 2376 /*25*/ MatZeroRows_SeqAIJ, 2377 MatLUFactorSymbolic_SeqAIJ, 2378 MatLUFactorNumeric_SeqAIJ, 2379 MatCholeskyFactorSymbolic_SeqAIJ, 2380 MatCholeskyFactorNumeric_SeqAIJ, 2381 /*30*/ MatSetUpPreallocation_SeqAIJ, 2382 MatILUFactorSymbolic_SeqAIJ, 2383 MatICCFactorSymbolic_SeqAIJ, 2384 MatGetArray_SeqAIJ, 2385 MatRestoreArray_SeqAIJ, 2386 /*35*/ MatDuplicate_SeqAIJ, 2387 0, 2388 0, 2389 MatILUFactor_SeqAIJ, 2390 0, 2391 /*40*/ MatAXPY_SeqAIJ, 2392 MatGetSubMatrices_SeqAIJ, 2393 MatIncreaseOverlap_SeqAIJ, 2394 MatGetValues_SeqAIJ, 2395 MatCopy_SeqAIJ, 2396 /*45*/ MatGetRowMax_SeqAIJ, 2397 MatScale_SeqAIJ, 2398 0, 2399 MatDiagonalSet_SeqAIJ, 2400 MatILUDTFactor_SeqAIJ, 2401 /*50*/ MatSetBlockSize_SeqAIJ, 2402 MatGetRowIJ_SeqAIJ, 2403 MatRestoreRowIJ_SeqAIJ, 2404 MatGetColumnIJ_SeqAIJ, 2405 MatRestoreColumnIJ_SeqAIJ, 2406 /*55*/ MatFDColoringCreate_SeqAIJ, 2407 0, 2408 0, 2409 MatPermute_SeqAIJ, 2410 0, 2411 /*60*/ 0, 2412 MatDestroy_SeqAIJ, 2413 MatView_SeqAIJ, 2414 0, 2415 0, 2416 /*65*/ 0, 2417 0, 2418 0, 2419 0, 2420 0, 2421 /*70*/ MatGetRowMaxAbs_SeqAIJ, 2422 0, 2423 MatSetColoring_SeqAIJ, 2424 #if defined(PETSC_HAVE_ADIC) 2425 MatSetValuesAdic_SeqAIJ, 2426 #else 2427 0, 2428 #endif 2429 MatSetValuesAdifor_SeqAIJ, 2430 /*75*/ 0, 2431 0, 2432 0, 2433 0, 2434 0, 2435 /*80*/ 0, 2436 0, 2437 0, 2438 0, 2439 MatLoad_SeqAIJ, 2440 /*85*/ MatIsSymmetric_SeqAIJ, 2441 MatIsHermitian_SeqAIJ, 2442 0, 2443 0, 2444 0, 2445 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2446 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2447 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2448 MatPtAP_Basic, 2449 MatPtAPSymbolic_SeqAIJ, 2450 /*95*/ MatPtAPNumeric_SeqAIJ, 2451 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2452 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2453 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2454 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2455 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2456 0, 2457 0, 2458 MatConjugate_SeqAIJ, 2459 0, 2460 /*105*/MatSetValuesRow_SeqAIJ, 2461 MatRealPart_SeqAIJ, 2462 MatImaginaryPart_SeqAIJ, 2463 0, 2464 0, 2465 /*110*/MatMatSolve_SeqAIJ, 2466 0, 2467 MatGetRowMin_SeqAIJ, 2468 0, 2469 MatMissingDiagonal_SeqAIJ 2470 }; 2471 2472 EXTERN_C_BEGIN 2473 #undef __FUNCT__ 2474 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2475 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2476 { 2477 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2478 PetscInt i,nz,n; 2479 2480 PetscFunctionBegin; 2481 2482 nz = aij->maxnz; 2483 n = mat->rmap.n; 2484 for (i=0; i<nz; i++) { 2485 aij->j[i] = indices[i]; 2486 } 2487 aij->nz = nz; 2488 for (i=0; i<n; i++) { 2489 aij->ilen[i] = aij->imax[i]; 2490 } 2491 2492 PetscFunctionReturn(0); 2493 } 2494 EXTERN_C_END 2495 2496 #undef __FUNCT__ 2497 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2498 /*@ 2499 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2500 in the matrix. 2501 2502 Input Parameters: 2503 + mat - the SeqAIJ matrix 2504 - indices - the column indices 2505 2506 Level: advanced 2507 2508 Notes: 2509 This can be called if you have precomputed the nonzero structure of the 2510 matrix and want to provide it to the matrix object to improve the performance 2511 of the MatSetValues() operation. 2512 2513 You MUST have set the correct numbers of nonzeros per row in the call to 2514 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2515 2516 MUST be called before any calls to MatSetValues(); 2517 2518 The indices should start with zero, not one. 2519 2520 @*/ 2521 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2522 { 2523 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2524 2525 PetscFunctionBegin; 2526 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2527 PetscValidPointer(indices,2); 2528 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2529 if (f) { 2530 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2531 } else { 2532 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2533 } 2534 PetscFunctionReturn(0); 2535 } 2536 2537 /* ----------------------------------------------------------------------------------------*/ 2538 2539 EXTERN_C_BEGIN 2540 #undef __FUNCT__ 2541 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2542 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2543 { 2544 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2545 PetscErrorCode ierr; 2546 size_t nz = aij->i[mat->rmap.n]; 2547 2548 PetscFunctionBegin; 2549 if (aij->nonew != 1) { 2550 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2551 } 2552 2553 /* allocate space for values if not already there */ 2554 if (!aij->saved_values) { 2555 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2556 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2557 } 2558 2559 /* copy values over */ 2560 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2561 PetscFunctionReturn(0); 2562 } 2563 EXTERN_C_END 2564 2565 #undef __FUNCT__ 2566 #define __FUNCT__ "MatStoreValues" 2567 /*@ 2568 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2569 example, reuse of the linear part of a Jacobian, while recomputing the 2570 nonlinear portion. 2571 2572 Collect on Mat 2573 2574 Input Parameters: 2575 . mat - the matrix (currently only AIJ matrices support this option) 2576 2577 Level: advanced 2578 2579 Common Usage, with SNESSolve(): 2580 $ Create Jacobian matrix 2581 $ Set linear terms into matrix 2582 $ Apply boundary conditions to matrix, at this time matrix must have 2583 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2584 $ boundary conditions again will not change the nonzero structure 2585 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2586 $ ierr = MatStoreValues(mat); 2587 $ Call SNESSetJacobian() with matrix 2588 $ In your Jacobian routine 2589 $ ierr = MatRetrieveValues(mat); 2590 $ Set nonlinear terms in matrix 2591 2592 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2593 $ // build linear portion of Jacobian 2594 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2595 $ ierr = MatStoreValues(mat); 2596 $ loop over nonlinear iterations 2597 $ ierr = MatRetrieveValues(mat); 2598 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2599 $ // call MatAssemblyBegin/End() on matrix 2600 $ Solve linear system with Jacobian 2601 $ endloop 2602 2603 Notes: 2604 Matrix must already be assemblied before calling this routine 2605 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2606 calling this routine. 2607 2608 When this is called multiple times it overwrites the previous set of stored values 2609 and does not allocated additional space. 2610 2611 .seealso: MatRetrieveValues() 2612 2613 @*/ 2614 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2615 { 2616 PetscErrorCode ierr,(*f)(Mat); 2617 2618 PetscFunctionBegin; 2619 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2620 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2621 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2622 2623 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2624 if (f) { 2625 ierr = (*f)(mat);CHKERRQ(ierr); 2626 } else { 2627 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2628 } 2629 PetscFunctionReturn(0); 2630 } 2631 2632 EXTERN_C_BEGIN 2633 #undef __FUNCT__ 2634 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2635 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2636 { 2637 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2638 PetscErrorCode ierr; 2639 PetscInt nz = aij->i[mat->rmap.n]; 2640 2641 PetscFunctionBegin; 2642 if (aij->nonew != 1) { 2643 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2644 } 2645 if (!aij->saved_values) { 2646 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2647 } 2648 /* copy values over */ 2649 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2650 PetscFunctionReturn(0); 2651 } 2652 EXTERN_C_END 2653 2654 #undef __FUNCT__ 2655 #define __FUNCT__ "MatRetrieveValues" 2656 /*@ 2657 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2658 example, reuse of the linear part of a Jacobian, while recomputing the 2659 nonlinear portion. 2660 2661 Collect on Mat 2662 2663 Input Parameters: 2664 . mat - the matrix (currently on AIJ matrices support this option) 2665 2666 Level: advanced 2667 2668 .seealso: MatStoreValues() 2669 2670 @*/ 2671 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2672 { 2673 PetscErrorCode ierr,(*f)(Mat); 2674 2675 PetscFunctionBegin; 2676 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2677 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2678 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2679 2680 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2681 if (f) { 2682 ierr = (*f)(mat);CHKERRQ(ierr); 2683 } else { 2684 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2685 } 2686 PetscFunctionReturn(0); 2687 } 2688 2689 2690 /* --------------------------------------------------------------------------------*/ 2691 #undef __FUNCT__ 2692 #define __FUNCT__ "MatCreateSeqAIJ" 2693 /*@C 2694 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2695 (the default parallel PETSc format). For good matrix assembly performance 2696 the user should preallocate the matrix storage by setting the parameter nz 2697 (or the array nnz). By setting these parameters accurately, performance 2698 during matrix assembly can be increased by more than a factor of 50. 2699 2700 Collective on MPI_Comm 2701 2702 Input Parameters: 2703 + comm - MPI communicator, set to PETSC_COMM_SELF 2704 . m - number of rows 2705 . n - number of columns 2706 . nz - number of nonzeros per row (same for all rows) 2707 - nnz - array containing the number of nonzeros in the various rows 2708 (possibly different for each row) or PETSC_NULL 2709 2710 Output Parameter: 2711 . A - the matrix 2712 2713 Notes: 2714 If nnz is given then nz is ignored 2715 2716 The AIJ format (also called the Yale sparse matrix format or 2717 compressed row storage), is fully compatible with standard Fortran 77 2718 storage. That is, the stored row and column indices can begin at 2719 either one (as in Fortran) or zero. See the users' manual for details. 2720 2721 Specify the preallocated storage with either nz or nnz (not both). 2722 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2723 allocation. For large problems you MUST preallocate memory or you 2724 will get TERRIBLE performance, see the users' manual chapter on matrices. 2725 2726 By default, this format uses inodes (identical nodes) when possible, to 2727 improve numerical efficiency of matrix-vector products and solves. We 2728 search for consecutive rows with the same nonzero structure, thereby 2729 reusing matrix information to achieve increased efficiency. 2730 2731 Options Database Keys: 2732 + -mat_no_inode - Do not use inodes 2733 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2734 - -mat_aij_oneindex - Internally use indexing starting at 1 2735 rather than 0. Note that when calling MatSetValues(), 2736 the user still MUST index entries starting at 0! 2737 2738 Level: intermediate 2739 2740 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2741 2742 @*/ 2743 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2744 { 2745 PetscErrorCode ierr; 2746 2747 PetscFunctionBegin; 2748 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2749 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2750 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2751 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2752 PetscFunctionReturn(0); 2753 } 2754 2755 #undef __FUNCT__ 2756 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2757 /*@C 2758 MatSeqAIJSetPreallocation - For good matrix assembly performance 2759 the user should preallocate the matrix storage by setting the parameter nz 2760 (or the array nnz). By setting these parameters accurately, performance 2761 during matrix assembly can be increased by more than a factor of 50. 2762 2763 Collective on MPI_Comm 2764 2765 Input Parameters: 2766 + B - The matrix 2767 . nz - number of nonzeros per row (same for all rows) 2768 - nnz - array containing the number of nonzeros in the various rows 2769 (possibly different for each row) or PETSC_NULL 2770 2771 Notes: 2772 If nnz is given then nz is ignored 2773 2774 The AIJ format (also called the Yale sparse matrix format or 2775 compressed row storage), is fully compatible with standard Fortran 77 2776 storage. That is, the stored row and column indices can begin at 2777 either one (as in Fortran) or zero. See the users' manual for details. 2778 2779 Specify the preallocated storage with either nz or nnz (not both). 2780 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2781 allocation. For large problems you MUST preallocate memory or you 2782 will get TERRIBLE performance, see the users' manual chapter on matrices. 2783 2784 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2785 entries or columns indices 2786 2787 By default, this format uses inodes (identical nodes) when possible, to 2788 improve numerical efficiency of matrix-vector products and solves. We 2789 search for consecutive rows with the same nonzero structure, thereby 2790 reusing matrix information to achieve increased efficiency. 2791 2792 Options Database Keys: 2793 + -mat_no_inode - Do not use inodes 2794 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2795 - -mat_aij_oneindex - Internally use indexing starting at 1 2796 rather than 0. Note that when calling MatSetValues(), 2797 the user still MUST index entries starting at 0! 2798 2799 Level: intermediate 2800 2801 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2802 2803 @*/ 2804 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2805 { 2806 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2807 2808 PetscFunctionBegin; 2809 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2810 if (f) { 2811 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2812 } 2813 PetscFunctionReturn(0); 2814 } 2815 2816 EXTERN_C_BEGIN 2817 #undef __FUNCT__ 2818 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2819 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2820 { 2821 Mat_SeqAIJ *b; 2822 PetscTruth skipallocation = PETSC_FALSE; 2823 PetscErrorCode ierr; 2824 PetscInt i; 2825 2826 PetscFunctionBegin; 2827 2828 if (nz == MAT_SKIP_ALLOCATION) { 2829 skipallocation = PETSC_TRUE; 2830 nz = 0; 2831 } 2832 2833 B->rmap.bs = B->cmap.bs = 1; 2834 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 2835 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 2836 2837 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2838 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2839 if (nnz) { 2840 for (i=0; i<B->rmap.n; i++) { 2841 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2842 if (nnz[i] > B->cmap.n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap.n); 2843 } 2844 } 2845 2846 B->preallocated = PETSC_TRUE; 2847 b = (Mat_SeqAIJ*)B->data; 2848 2849 if (!skipallocation) { 2850 ierr = PetscMalloc2(B->rmap.n,PetscInt,&b->imax,B->rmap.n,PetscInt,&b->ilen);CHKERRQ(ierr); 2851 ierr = PetscLogObjectMemory(B,2*B->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 2852 if (!nnz) { 2853 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2854 else if (nz <= 0) nz = 1; 2855 for (i=0; i<B->rmap.n; i++) b->imax[i] = nz; 2856 nz = nz*B->rmap.n; 2857 } else { 2858 nz = 0; 2859 for (i=0; i<B->rmap.n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2860 } 2861 2862 /* b->ilen will count nonzeros in each row so far. */ 2863 for (i=0; i<B->rmap.n; i++) { b->ilen[i] = 0;} 2864 2865 /* allocate the matrix space */ 2866 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.n+1,PetscInt,&b->i);CHKERRQ(ierr); 2867 ierr = PetscLogObjectMemory(B,(B->rmap.n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 2868 b->i[0] = 0; 2869 for (i=1; i<B->rmap.n+1; i++) { 2870 b->i[i] = b->i[i-1] + b->imax[i-1]; 2871 } 2872 b->singlemalloc = PETSC_TRUE; 2873 b->free_a = PETSC_TRUE; 2874 b->free_ij = PETSC_TRUE; 2875 } else { 2876 b->free_a = PETSC_FALSE; 2877 b->free_ij = PETSC_FALSE; 2878 } 2879 2880 b->nz = 0; 2881 b->maxnz = nz; 2882 B->info.nz_unneeded = (double)b->maxnz; 2883 PetscFunctionReturn(0); 2884 } 2885 EXTERN_C_END 2886 2887 #undef __FUNCT__ 2888 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2889 /*@C 2890 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2891 2892 Input Parameters: 2893 + B - the matrix 2894 . i - the indices into j for the start of each row (starts with zero) 2895 . j - the column indices for each row (starts with zero) these must be sorted for each row 2896 - v - optional values in the matrix 2897 2898 Contributed by: Lisandro Dalchin 2899 2900 Level: developer 2901 2902 .keywords: matrix, aij, compressed row, sparse, sequential 2903 2904 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2905 @*/ 2906 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2907 { 2908 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2909 PetscErrorCode ierr; 2910 2911 PetscFunctionBegin; 2912 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2913 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2914 if (f) { 2915 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2916 } 2917 PetscFunctionReturn(0); 2918 } 2919 2920 EXTERN_C_BEGIN 2921 #undef __FUNCT__ 2922 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2923 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 2924 { 2925 PetscInt i; 2926 PetscInt m,n; 2927 PetscInt nz; 2928 PetscInt *nnz, nz_max = 0; 2929 PetscScalar *values; 2930 PetscErrorCode ierr; 2931 2932 PetscFunctionBegin; 2933 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2934 2935 if (Ii[0]) { 2936 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 2937 } 2938 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2939 for(i = 0; i < m; i++) { 2940 nz = Ii[i+1]- Ii[i]; 2941 nz_max = PetscMax(nz_max, nz); 2942 if (nz < 0) { 2943 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2944 } 2945 nnz[i] = nz; 2946 } 2947 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2948 ierr = PetscFree(nnz);CHKERRQ(ierr); 2949 2950 if (v) { 2951 values = (PetscScalar*) v; 2952 } else { 2953 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2954 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2955 } 2956 2957 for(i = 0; i < m; i++) { 2958 nz = Ii[i+1] - Ii[i]; 2959 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2960 } 2961 2962 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2963 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2964 2965 if (!v) { 2966 ierr = PetscFree(values);CHKERRQ(ierr); 2967 } 2968 PetscFunctionReturn(0); 2969 } 2970 EXTERN_C_END 2971 2972 /*MC 2973 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2974 based on compressed sparse row format. 2975 2976 Options Database Keys: 2977 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2978 2979 Level: beginner 2980 2981 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2982 M*/ 2983 2984 EXTERN_C_BEGIN 2985 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 2986 EXTERN_C_END 2987 2988 EXTERN_C_BEGIN 2989 #undef __FUNCT__ 2990 #define __FUNCT__ "MatCreate_SeqAIJ" 2991 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 2992 { 2993 Mat_SeqAIJ *b; 2994 PetscErrorCode ierr; 2995 PetscMPIInt size; 2996 2997 PetscFunctionBegin; 2998 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2999 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3000 3001 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3002 B->data = (void*)b; 3003 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3004 B->factor = 0; 3005 B->mapping = 0; 3006 b->row = 0; 3007 b->col = 0; 3008 b->icol = 0; 3009 b->reallocs = 0; 3010 b->ignorezeroentries = PETSC_FALSE; 3011 b->roworiented = PETSC_TRUE; 3012 b->nonew = 0; 3013 b->diag = 0; 3014 b->solve_work = 0; 3015 B->spptr = 0; 3016 b->saved_values = 0; 3017 b->idiag = 0; 3018 b->ssor = 0; 3019 b->keepzeroedrows = PETSC_FALSE; 3020 b->xtoy = 0; 3021 b->XtoY = 0; 3022 b->compressedrow.use = PETSC_FALSE; 3023 b->compressedrow.nrows = B->rmap.n; 3024 b->compressedrow.i = PETSC_NULL; 3025 b->compressedrow.rindex = PETSC_NULL; 3026 b->compressedrow.checked = PETSC_FALSE; 3027 B->same_nonzero = PETSC_FALSE; 3028 3029 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3030 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3031 "MatSeqAIJSetColumnIndices_SeqAIJ", 3032 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3033 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3034 "MatStoreValues_SeqAIJ", 3035 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3036 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3037 "MatRetrieveValues_SeqAIJ", 3038 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3039 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3040 "MatConvert_SeqAIJ_SeqSBAIJ", 3041 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3042 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3043 "MatConvert_SeqAIJ_SeqBAIJ", 3044 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3045 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3046 "MatConvert_SeqAIJ_SeqCSRPERM", 3047 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3048 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3049 "MatConvert_SeqAIJ_SeqCRL", 3050 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3051 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3052 "MatIsTranspose_SeqAIJ", 3053 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3054 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3055 "MatIsHermitianTranspose_SeqAIJ", 3056 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3057 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3058 "MatSeqAIJSetPreallocation_SeqAIJ", 3059 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3060 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3061 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3062 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3063 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3064 "MatReorderForNonzeroDiagonal_SeqAIJ", 3065 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3066 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 3067 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3068 PetscFunctionReturn(0); 3069 } 3070 EXTERN_C_END 3071 3072 #undef __FUNCT__ 3073 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3074 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3075 { 3076 Mat C; 3077 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3078 PetscErrorCode ierr; 3079 PetscInt i,m = A->rmap.n; 3080 3081 PetscFunctionBegin; 3082 *B = 0; 3083 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 3084 ierr = MatSetSizes(C,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 3085 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 3086 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3087 3088 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 3089 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 3090 3091 c = (Mat_SeqAIJ*)C->data; 3092 3093 C->factor = A->factor; 3094 3095 c->row = 0; 3096 c->col = 0; 3097 c->icol = 0; 3098 c->reallocs = 0; 3099 3100 C->assembled = PETSC_TRUE; 3101 3102 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3103 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3104 for (i=0; i<m; i++) { 3105 c->imax[i] = a->imax[i]; 3106 c->ilen[i] = a->ilen[i]; 3107 } 3108 3109 /* allocate the matrix space */ 3110 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3111 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3112 c->singlemalloc = PETSC_TRUE; 3113 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3114 if (m > 0) { 3115 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3116 if (cpvalues == MAT_COPY_VALUES) { 3117 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3118 } else { 3119 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3120 } 3121 } 3122 3123 c->ignorezeroentries = a->ignorezeroentries; 3124 c->roworiented = a->roworiented; 3125 c->nonew = a->nonew; 3126 if (a->diag) { 3127 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3128 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3129 for (i=0; i<m; i++) { 3130 c->diag[i] = a->diag[i]; 3131 } 3132 } else c->diag = 0; 3133 c->solve_work = 0; 3134 c->saved_values = 0; 3135 c->idiag = 0; 3136 c->ssor = 0; 3137 c->keepzeroedrows = a->keepzeroedrows; 3138 c->free_a = PETSC_TRUE; 3139 c->free_ij = PETSC_TRUE; 3140 c->xtoy = 0; 3141 c->XtoY = 0; 3142 3143 c->nz = a->nz; 3144 c->maxnz = a->maxnz; 3145 C->preallocated = PETSC_TRUE; 3146 3147 c->compressedrow.use = a->compressedrow.use; 3148 c->compressedrow.nrows = a->compressedrow.nrows; 3149 c->compressedrow.checked = a->compressedrow.checked; 3150 if ( a->compressedrow.checked && a->compressedrow.use){ 3151 i = a->compressedrow.nrows; 3152 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3153 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3154 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3155 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3156 } else { 3157 c->compressedrow.use = PETSC_FALSE; 3158 c->compressedrow.i = PETSC_NULL; 3159 c->compressedrow.rindex = PETSC_NULL; 3160 } 3161 C->same_nonzero = A->same_nonzero; 3162 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3163 3164 *B = C; 3165 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 3166 PetscFunctionReturn(0); 3167 } 3168 3169 #undef __FUNCT__ 3170 #define __FUNCT__ "MatLoad_SeqAIJ" 3171 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 3172 { 3173 Mat_SeqAIJ *a; 3174 Mat B; 3175 PetscErrorCode ierr; 3176 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3177 int fd; 3178 PetscMPIInt size; 3179 MPI_Comm comm; 3180 3181 PetscFunctionBegin; 3182 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3183 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3184 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3185 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3186 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3187 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3188 M = header[1]; N = header[2]; nz = header[3]; 3189 3190 if (nz < 0) { 3191 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3192 } 3193 3194 /* read in row lengths */ 3195 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3196 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3197 3198 /* check if sum of rowlengths is same as nz */ 3199 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3200 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3201 3202 /* create our matrix */ 3203 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3204 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3205 ierr = MatSetType(B,type);CHKERRQ(ierr); 3206 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3207 a = (Mat_SeqAIJ*)B->data; 3208 3209 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3210 3211 /* read in nonzero values */ 3212 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3213 3214 /* set matrix "i" values */ 3215 a->i[0] = 0; 3216 for (i=1; i<= M; i++) { 3217 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3218 a->ilen[i-1] = rowlengths[i-1]; 3219 } 3220 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3221 3222 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3223 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3224 *A = B; 3225 PetscFunctionReturn(0); 3226 } 3227 3228 #undef __FUNCT__ 3229 #define __FUNCT__ "MatEqual_SeqAIJ" 3230 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3231 { 3232 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3233 PetscErrorCode ierr; 3234 3235 PetscFunctionBegin; 3236 /* If the matrix dimensions are not equal,or no of nonzeros */ 3237 if ((A->rmap.n != B->rmap.n) || (A->cmap.n != B->cmap.n) ||(a->nz != b->nz)) { 3238 *flg = PETSC_FALSE; 3239 PetscFunctionReturn(0); 3240 } 3241 3242 /* if the a->i are the same */ 3243 ierr = PetscMemcmp(a->i,b->i,(A->rmap.n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3244 if (!*flg) PetscFunctionReturn(0); 3245 3246 /* if a->j are the same */ 3247 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3248 if (!*flg) PetscFunctionReturn(0); 3249 3250 /* if a->a are the same */ 3251 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3252 3253 PetscFunctionReturn(0); 3254 3255 } 3256 3257 #undef __FUNCT__ 3258 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3259 /*@ 3260 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3261 provided by the user. 3262 3263 Collective on MPI_Comm 3264 3265 Input Parameters: 3266 + comm - must be an MPI communicator of size 1 3267 . m - number of rows 3268 . n - number of columns 3269 . i - row indices 3270 . j - column indices 3271 - a - matrix values 3272 3273 Output Parameter: 3274 . mat - the matrix 3275 3276 Level: intermediate 3277 3278 Notes: 3279 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3280 once the matrix is destroyed 3281 3282 You cannot set new nonzero locations into this matrix, that will generate an error. 3283 3284 The i and j indices are 0 based 3285 3286 The format which is used for the sparse matrix input, is equivalent to a 3287 row-major ordering.. i.e for the following matrix, the input data expected is 3288 as shown: 3289 3290 1 0 0 3291 2 0 3 3292 4 5 6 3293 3294 i = {0,1,3,6} [size = nrow+1 = 3+1] 3295 j = {0,0,2,0,1,2} [size = nz = 6] 3296 v = {1,2,3,4,5,6} [size = nz = 6] 3297 3298 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3299 3300 @*/ 3301 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3302 { 3303 PetscErrorCode ierr; 3304 PetscInt ii; 3305 Mat_SeqAIJ *aij; 3306 3307 PetscFunctionBegin; 3308 if (i[0]) { 3309 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3310 } 3311 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3312 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3313 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3314 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3315 aij = (Mat_SeqAIJ*)(*mat)->data; 3316 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3317 3318 aij->i = i; 3319 aij->j = j; 3320 aij->a = a; 3321 aij->singlemalloc = PETSC_FALSE; 3322 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3323 aij->free_a = PETSC_FALSE; 3324 aij->free_ij = PETSC_FALSE; 3325 3326 for (ii=0; ii<m; ii++) { 3327 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3328 #if defined(PETSC_USE_DEBUG) 3329 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3330 #endif 3331 } 3332 #if defined(PETSC_USE_DEBUG) 3333 for (ii=0; ii<aij->i[m]; ii++) { 3334 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3335 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3336 } 3337 #endif 3338 3339 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3340 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3341 PetscFunctionReturn(0); 3342 } 3343 3344 #undef __FUNCT__ 3345 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3346 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3347 { 3348 PetscErrorCode ierr; 3349 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3350 3351 PetscFunctionBegin; 3352 if (coloring->ctype == IS_COLORING_GLOBAL) { 3353 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3354 a->coloring = coloring; 3355 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3356 PetscInt i,*larray; 3357 ISColoring ocoloring; 3358 ISColoringValue *colors; 3359 3360 /* set coloring for diagonal portion */ 3361 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3362 for (i=0; i<A->cmap.n; i++) { 3363 larray[i] = i; 3364 } 3365 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap.n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3366 ierr = PetscMalloc((A->cmap.n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3367 for (i=0; i<A->cmap.n; i++) { 3368 colors[i] = coloring->colors[larray[i]]; 3369 } 3370 ierr = PetscFree(larray);CHKERRQ(ierr); 3371 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap.n,colors,&ocoloring);CHKERRQ(ierr); 3372 a->coloring = ocoloring; 3373 } 3374 PetscFunctionReturn(0); 3375 } 3376 3377 #if defined(PETSC_HAVE_ADIC) 3378 EXTERN_C_BEGIN 3379 #include "adic/ad_utils.h" 3380 EXTERN_C_END 3381 3382 #undef __FUNCT__ 3383 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3384 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3385 { 3386 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3387 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3388 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3389 ISColoringValue *color; 3390 3391 PetscFunctionBegin; 3392 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3393 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3394 color = a->coloring->colors; 3395 /* loop over rows */ 3396 for (i=0; i<m; i++) { 3397 nz = ii[i+1] - ii[i]; 3398 /* loop over columns putting computed value into matrix */ 3399 for (j=0; j<nz; j++) { 3400 *v++ = values[color[*jj++]]; 3401 } 3402 values += nlen; /* jump to next row of derivatives */ 3403 } 3404 PetscFunctionReturn(0); 3405 } 3406 #endif 3407 3408 #undef __FUNCT__ 3409 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3410 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3411 { 3412 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3413 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j; 3414 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3415 ISColoringValue *color; 3416 3417 PetscFunctionBegin; 3418 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3419 color = a->coloring->colors; 3420 /* loop over rows */ 3421 for (i=0; i<m; i++) { 3422 nz = ii[i+1] - ii[i]; 3423 /* loop over columns putting computed value into matrix */ 3424 for (j=0; j<nz; j++) { 3425 *v++ = values[color[*jj++]]; 3426 } 3427 values += nl; /* jump to next row of derivatives */ 3428 } 3429 PetscFunctionReturn(0); 3430 } 3431 3432 /* 3433 Special version for direct calls from Fortran 3434 */ 3435 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3436 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3437 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3438 #define matsetvaluesseqaij_ matsetvaluesseqaij 3439 #endif 3440 3441 /* Change these macros so can be used in void function */ 3442 #undef CHKERRQ 3443 #define CHKERRQ(ierr) CHKERRABORT(A->comm,ierr) 3444 #undef SETERRQ2 3445 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(A->comm,ierr) 3446 3447 EXTERN_C_BEGIN 3448 #undef __FUNCT__ 3449 #define __FUNCT__ "matsetvaluesseqaij_" 3450 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3451 { 3452 Mat A = *AA; 3453 PetscInt m = *mm, n = *nn; 3454 InsertMode is = *isis; 3455 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3456 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3457 PetscInt *imax,*ai,*ailen; 3458 PetscErrorCode ierr; 3459 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3460 PetscScalar *ap,value,*aa; 3461 PetscTruth ignorezeroentries = a->ignorezeroentries; 3462 PetscTruth roworiented = a->roworiented; 3463 3464 PetscFunctionBegin; 3465 MatPreallocated(A); 3466 imax = a->imax; 3467 ai = a->i; 3468 ailen = a->ilen; 3469 aj = a->j; 3470 aa = a->a; 3471 3472 for (k=0; k<m; k++) { /* loop over added rows */ 3473 row = im[k]; 3474 if (row < 0) continue; 3475 #if defined(PETSC_USE_DEBUG) 3476 if (row >= A->rmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3477 #endif 3478 rp = aj + ai[row]; ap = aa + ai[row]; 3479 rmax = imax[row]; nrow = ailen[row]; 3480 low = 0; 3481 high = nrow; 3482 for (l=0; l<n; l++) { /* loop over added columns */ 3483 if (in[l] < 0) continue; 3484 #if defined(PETSC_USE_DEBUG) 3485 if (in[l] >= A->cmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3486 #endif 3487 col = in[l]; 3488 if (roworiented) { 3489 value = v[l + k*n]; 3490 } else { 3491 value = v[k + l*m]; 3492 } 3493 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3494 3495 if (col <= lastcol) low = 0; else high = nrow; 3496 lastcol = col; 3497 while (high-low > 5) { 3498 t = (low+high)/2; 3499 if (rp[t] > col) high = t; 3500 else low = t; 3501 } 3502 for (i=low; i<high; i++) { 3503 if (rp[i] > col) break; 3504 if (rp[i] == col) { 3505 if (is == ADD_VALUES) ap[i] += value; 3506 else ap[i] = value; 3507 goto noinsert; 3508 } 3509 } 3510 if (value == 0.0 && ignorezeroentries) goto noinsert; 3511 if (nonew == 1) goto noinsert; 3512 if (nonew == -1) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3513 MatSeqXAIJReallocateAIJ(A,A->rmap.n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3514 N = nrow++ - 1; a->nz++; high++; 3515 /* shift up all the later entries in this row */ 3516 for (ii=N; ii>=i; ii--) { 3517 rp[ii+1] = rp[ii]; 3518 ap[ii+1] = ap[ii]; 3519 } 3520 rp[i] = col; 3521 ap[i] = value; 3522 noinsert:; 3523 low = i + 1; 3524 } 3525 ailen[row] = nrow; 3526 } 3527 A->same_nonzero = PETSC_FALSE; 3528 PetscFunctionReturnVoid(); 3529 } 3530 EXTERN_C_END 3531