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