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