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