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