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