1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the ADJ adjacency list matrix data-structure. 5 */ 6 #include "../src/mat/impls/adj/mpi/mpiadj.h" 7 #include "petscsys.h" 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatView_MPIAdj_ASCII" 11 PetscErrorCode MatView_MPIAdj_ASCII(Mat A,PetscViewer viewer) 12 { 13 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 14 PetscErrorCode ierr; 15 PetscInt i,j,m = A->rmap->n; 16 const char *name; 17 PetscViewerFormat format; 18 19 PetscFunctionBegin; 20 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 21 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 22 if (format == PETSC_VIEWER_ASCII_INFO) { 23 PetscFunctionReturn(0); 24 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 25 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 26 } else { 27 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 28 for (i=0; i<m; i++) { 29 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"row %D:",i+A->rmap->rstart);CHKERRQ(ierr); 30 for (j=a->i[i]; j<a->i[i+1]; j++) { 31 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %D ",a->j[j]);CHKERRQ(ierr); 32 } 33 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 34 } 35 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 36 } 37 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 #undef __FUNCT__ 42 #define __FUNCT__ "MatView_MPIAdj" 43 PetscErrorCode MatView_MPIAdj(Mat A,PetscViewer viewer) 44 { 45 PetscErrorCode ierr; 46 PetscTruth iascii; 47 48 PetscFunctionBegin; 49 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 50 if (iascii) { 51 ierr = MatView_MPIAdj_ASCII(A,viewer);CHKERRQ(ierr); 52 } else { 53 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by MPIAdj",((PetscObject)viewer)->type_name); 54 } 55 PetscFunctionReturn(0); 56 } 57 58 #undef __FUNCT__ 59 #define __FUNCT__ "MatDestroy_MPIAdj" 60 PetscErrorCode MatDestroy_MPIAdj(Mat mat) 61 { 62 Mat_MPIAdj *a = (Mat_MPIAdj*)mat->data; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 #if defined(PETSC_USE_LOG) 67 PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D, NZ=%D",mat->rmap->n,mat->cmap->n,a->nz); 68 #endif 69 ierr = PetscFree(a->diag);CHKERRQ(ierr); 70 if (a->freeaij) { 71 ierr = PetscFree(a->i);CHKERRQ(ierr); 72 ierr = PetscFree(a->j);CHKERRQ(ierr); 73 ierr = PetscFree(a->values);CHKERRQ(ierr); 74 } else { 75 if (a->i) free(a->i); 76 if (a->j) free(a->j); 77 } 78 ierr = PetscFree(a);CHKERRQ(ierr); 79 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 80 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMPIAdjSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "MatSetOption_MPIAdj" 86 PetscErrorCode MatSetOption_MPIAdj(Mat A,MatOption op,PetscTruth flg) 87 { 88 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 switch (op) { 93 case MAT_SYMMETRIC: 94 case MAT_STRUCTURALLY_SYMMETRIC: 95 case MAT_HERMITIAN: 96 a->symmetric = flg; 97 break; 98 case MAT_SYMMETRY_ETERNAL: 99 break; 100 default: 101 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 102 break; 103 } 104 PetscFunctionReturn(0); 105 } 106 107 108 /* 109 Adds diagonal pointers to sparse matrix structure. 110 */ 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "MatMarkDiagonal_MPIAdj" 114 PetscErrorCode MatMarkDiagonal_MPIAdj(Mat A) 115 { 116 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 117 PetscErrorCode ierr; 118 PetscInt i,j,m = A->rmap->n; 119 120 PetscFunctionBegin; 121 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 122 ierr = PetscLogObjectMemory(A,m*sizeof(PetscInt));CHKERRQ(ierr); 123 for (i=0; i<A->rmap->n; i++) { 124 for (j=a->i[i]; j<a->i[i+1]; j++) { 125 if (a->j[j] == i) { 126 a->diag[i] = j; 127 break; 128 } 129 } 130 } 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "MatGetRow_MPIAdj" 136 PetscErrorCode MatGetRow_MPIAdj(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 137 { 138 Mat_MPIAdj *a = (Mat_MPIAdj*)A->data; 139 PetscInt *itmp; 140 141 PetscFunctionBegin; 142 row -= A->rmap->rstart; 143 144 if (row < 0 || row >= A->rmap->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 145 146 *nz = a->i[row+1] - a->i[row]; 147 if (v) *v = PETSC_NULL; 148 if (idx) { 149 itmp = a->j + a->i[row]; 150 if (*nz) { 151 *idx = itmp; 152 } 153 else *idx = 0; 154 } 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "MatRestoreRow_MPIAdj" 160 PetscErrorCode MatRestoreRow_MPIAdj(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 161 { 162 PetscFunctionBegin; 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "MatEqual_MPIAdj" 168 PetscErrorCode MatEqual_MPIAdj(Mat A,Mat B,PetscTruth* flg) 169 { 170 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data,*b = (Mat_MPIAdj *)B->data; 171 PetscErrorCode ierr; 172 PetscTruth flag; 173 174 PetscFunctionBegin; 175 /* If the matrix dimensions are not equal,or no of nonzeros */ 176 if ((A->rmap->n != B->rmap->n) ||(a->nz != b->nz)) { 177 flag = PETSC_FALSE; 178 } 179 180 /* if the a->i are the same */ 181 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),&flag);CHKERRQ(ierr); 182 183 /* if a->j are the same */ 184 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),&flag);CHKERRQ(ierr); 185 186 ierr = MPI_Allreduce(&flag,flg,1,MPI_INT,MPI_LAND,((PetscObject)A)->comm);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 #undef __FUNCT__ 191 #define __FUNCT__ "MatGetRowIJ_MPIAdj" 192 PetscErrorCode MatGetRowIJ_MPIAdj(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 193 { 194 PetscErrorCode ierr; 195 PetscMPIInt size; 196 PetscInt i; 197 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 198 199 PetscFunctionBegin; 200 ierr = MPI_Comm_size(((PetscObject)A)->comm,&size);CHKERRQ(ierr); 201 if (size > 1) {*done = PETSC_FALSE; PetscFunctionReturn(0);} 202 *m = A->rmap->n; 203 *ia = a->i; 204 *ja = a->j; 205 *done = PETSC_TRUE; 206 if (oshift) { 207 for (i=0; i<(*ia)[*m]; i++) { 208 (*ja)[i]++; 209 } 210 for (i=0; i<=(*m); i++) (*ia)[i]++; 211 } 212 PetscFunctionReturn(0); 213 } 214 215 #undef __FUNCT__ 216 #define __FUNCT__ "MatRestoreRowIJ_MPIAdj" 217 PetscErrorCode MatRestoreRowIJ_MPIAdj(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 218 { 219 PetscInt i; 220 Mat_MPIAdj *a = (Mat_MPIAdj *)A->data; 221 222 PetscFunctionBegin; 223 if (ia && a->i != *ia) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ia passed back is not one obtained with MatGetRowIJ()"); 224 if (ja && a->j != *ja) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ja passed back is not one obtained with MatGetRowIJ()"); 225 if (oshift) { 226 for (i=0; i<=(*m); i++) (*ia)[i]--; 227 for (i=0; i<(*ia)[*m]; i++) { 228 (*ja)[i]--; 229 } 230 } 231 PetscFunctionReturn(0); 232 } 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "MatConvertFrom_MPIAdj" 236 PetscErrorCode PETSCMAT_DLLEXPORT MatConvertFrom_MPIAdj(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 237 { 238 Mat B; 239 PetscErrorCode ierr; 240 PetscInt i,m,N,nzeros = 0,*ia,*ja,len,rstart,cnt,j,*a; 241 const PetscInt *rj; 242 const PetscScalar *ra; 243 MPI_Comm comm; 244 245 PetscFunctionBegin; 246 ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 247 ierr = MatGetLocalSize(A,&m,PETSC_NULL);CHKERRQ(ierr); 248 ierr = MatGetOwnershipRange(A,&rstart,PETSC_NULL);CHKERRQ(ierr); 249 250 /* count the number of nonzeros per row */ 251 for (i=0; i<m; i++) { 252 ierr = MatGetRow(A,i+rstart,&len,&rj,PETSC_NULL);CHKERRQ(ierr); 253 for (j=0; j<len; j++) { 254 if (rj[j] == i+rstart) {len--; break;} /* don't count diagonal */ 255 } 256 ierr = MatRestoreRow(A,i+rstart,&len,&rj,PETSC_NULL);CHKERRQ(ierr); 257 nzeros += len; 258 } 259 260 /* malloc space for nonzeros */ 261 ierr = PetscMalloc((nzeros+1)*sizeof(PetscInt),&a);CHKERRQ(ierr); 262 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&ia);CHKERRQ(ierr); 263 ierr = PetscMalloc((nzeros+1)*sizeof(PetscInt),&ja);CHKERRQ(ierr); 264 265 nzeros = 0; 266 ia[0] = 0; 267 for (i=0; i<m; i++) { 268 ierr = MatGetRow(A,i+rstart,&len,&rj,&ra);CHKERRQ(ierr); 269 cnt = 0; 270 for (j=0; j<len; j++) { 271 if (rj[j] != i+rstart) { /* if not diagonal */ 272 a[nzeros+cnt] = (PetscInt) PetscAbsScalar(ra[j]); 273 ja[nzeros+cnt++] = rj[j]; 274 } 275 } 276 ierr = MatRestoreRow(A,i+rstart,&len,&rj,&ra);CHKERRQ(ierr); 277 nzeros += cnt; 278 ia[i+1] = nzeros; 279 } 280 281 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 282 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 283 ierr = MatSetSizes(B,m,PETSC_DETERMINE,PETSC_DETERMINE,N);CHKERRQ(ierr); 284 ierr = MatSetType(B,type);CHKERRQ(ierr); 285 ierr = MatMPIAdjSetPreallocation(B,ia,ja,a);CHKERRQ(ierr); 286 287 if (reuse == MAT_REUSE_MATRIX) { 288 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 289 } else { 290 *newmat = B; 291 } 292 PetscFunctionReturn(0); 293 } 294 295 /* -------------------------------------------------------------------*/ 296 static struct _MatOps MatOps_Values = {0, 297 MatGetRow_MPIAdj, 298 MatRestoreRow_MPIAdj, 299 0, 300 /* 4*/ 0, 301 0, 302 0, 303 0, 304 0, 305 0, 306 /*10*/ 0, 307 0, 308 0, 309 0, 310 0, 311 /*15*/ 0, 312 MatEqual_MPIAdj, 313 0, 314 0, 315 0, 316 /*20*/ 0, 317 0, 318 MatSetOption_MPIAdj, 319 0, 320 /*24*/ 0, 321 0, 322 0, 323 0, 324 0, 325 /*29*/ 0, 326 0, 327 0, 328 0, 329 0, 330 /*34*/ 0, 331 0, 332 0, 333 0, 334 0, 335 /*39*/ 0, 336 0, 337 0, 338 0, 339 0, 340 /*44*/ 0, 341 0, 342 0, 343 0, 344 0, 345 /*49*/ 0, 346 MatGetRowIJ_MPIAdj, 347 MatRestoreRowIJ_MPIAdj, 348 0, 349 0, 350 /*54*/ 0, 351 0, 352 0, 353 0, 354 0, 355 /*59*/ 0, 356 MatDestroy_MPIAdj, 357 MatView_MPIAdj, 358 MatConvertFrom_MPIAdj, 359 0, 360 /*64*/ 0, 361 0, 362 0, 363 0, 364 0, 365 /*69*/ 0, 366 0, 367 0, 368 0, 369 0, 370 /*74*/ 0, 371 0, 372 0, 373 0, 374 0, 375 /*79*/ 0, 376 0, 377 0, 378 0, 379 0, 380 /*84*/ 0, 381 0, 382 0, 383 0, 384 0, 385 /*89*/ 0, 386 0, 387 0, 388 0, 389 0, 390 /*94*/ 0, 391 0, 392 0, 393 0}; 394 395 EXTERN_C_BEGIN 396 #undef __FUNCT__ 397 #define __FUNCT__ "MatMPIAdjSetPreallocation_MPIAdj" 398 PetscErrorCode PETSCMAT_DLLEXPORT MatMPIAdjSetPreallocation_MPIAdj(Mat B,PetscInt *i,PetscInt *j,PetscInt *values) 399 { 400 Mat_MPIAdj *b = (Mat_MPIAdj *)B->data; 401 PetscErrorCode ierr; 402 #if defined(PETSC_USE_DEBUG) 403 PetscInt ii; 404 #endif 405 406 PetscFunctionBegin; 407 ierr = PetscMapSetBlockSize(B->rmap,1);CHKERRQ(ierr); 408 ierr = PetscMapSetBlockSize(B->cmap,1);CHKERRQ(ierr); 409 ierr = PetscMapSetUp(B->rmap);CHKERRQ(ierr); 410 ierr = PetscMapSetUp(B->cmap);CHKERRQ(ierr); 411 412 #if defined(PETSC_USE_DEBUG) 413 if (i[0] != 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"First i[] index must be zero, instead it is %D\n",i[0]); 414 for (ii=1; ii<B->rmap->n; ii++) { 415 if (i[ii] < 0 || i[ii] < i[ii-1]) { 416 SETERRQ4(PETSC_ERR_ARG_OUTOFRANGE,"i[%D]=%D index is out of range: i[%D]=%D",ii,i[ii],ii-1,i[ii-1]); 417 } 418 } 419 for (ii=0; ii<i[B->rmap->n]; ii++) { 420 if (j[ii] < 0 || j[ii] >= B->cmap->N) { 421 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index %D out of range %D\n",ii,j[ii]); 422 } 423 } 424 #endif 425 B->preallocated = PETSC_TRUE; 426 427 b->j = j; 428 b->i = i; 429 b->values = values; 430 431 b->nz = i[B->rmap->n]; 432 b->diag = 0; 433 b->symmetric = PETSC_FALSE; 434 b->freeaij = PETSC_TRUE; 435 436 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 437 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 EXTERN_C_END 441 442 /*MC 443 MATMPIADJ - MATMPIADJ = "mpiadj" - A matrix type to be used for distributed adjacency matrices, 444 intended for use constructing orderings and partitionings. 445 446 Level: beginner 447 448 .seealso: MatCreateMPIAdj 449 M*/ 450 451 EXTERN_C_BEGIN 452 #undef __FUNCT__ 453 #define __FUNCT__ "MatCreate_MPIAdj" 454 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_MPIAdj(Mat B) 455 { 456 Mat_MPIAdj *b; 457 PetscErrorCode ierr; 458 PetscMPIInt size,rank; 459 460 PetscFunctionBegin; 461 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 462 ierr = MPI_Comm_rank(((PetscObject)B)->comm,&rank);CHKERRQ(ierr); 463 464 ierr = PetscNewLog(B,Mat_MPIAdj,&b);CHKERRQ(ierr); 465 B->data = (void*)b; 466 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 467 B->mapping = 0; 468 B->assembled = PETSC_FALSE; 469 470 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMPIAdjSetPreallocation_C", 471 "MatMPIAdjSetPreallocation_MPIAdj", 472 MatMPIAdjSetPreallocation_MPIAdj);CHKERRQ(ierr); 473 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMPIADJ);CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476 EXTERN_C_END 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "MatMPIAdjSetPreallocation" 480 /*@C 481 MatMPIAdjSetPreallocation - Sets the array used for storing the matrix elements 482 483 Collective on MPI_Comm 484 485 Input Parameters: 486 + A - the matrix 487 . i - the indices into j for the start of each row 488 . j - the column indices for each row (sorted for each row). 489 The indices in i and j start with zero (NOT with one). 490 - values - [optional] edge weights 491 492 Level: intermediate 493 494 .seealso: MatCreate(), MatCreateMPIAdj(), MatSetValues() 495 @*/ 496 PetscErrorCode PETSCMAT_DLLEXPORT MatMPIAdjSetPreallocation(Mat B,PetscInt *i,PetscInt *j,PetscInt *values) 497 { 498 PetscErrorCode ierr,(*f)(Mat,PetscInt*,PetscInt*,PetscInt*); 499 500 PetscFunctionBegin; 501 ierr = PetscObjectQueryFunction((PetscObject)B,"MatMPIAdjSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 502 if (f) { 503 ierr = (*f)(B,i,j,values);CHKERRQ(ierr); 504 } 505 PetscFunctionReturn(0); 506 } 507 508 #undef __FUNCT__ 509 #define __FUNCT__ "MatCreateMPIAdj" 510 /*@C 511 MatCreateMPIAdj - Creates a sparse matrix representing an adjacency list. 512 The matrix does not have numerical values associated with it, but is 513 intended for ordering (to reduce bandwidth etc) and partitioning. 514 515 Collective on MPI_Comm 516 517 Input Parameters: 518 + comm - MPI communicator 519 . m - number of local rows 520 . N - number of global columns 521 . i - the indices into j for the start of each row 522 . j - the column indices for each row (sorted for each row). 523 The indices in i and j start with zero (NOT with one). 524 - values -[optional] edge weights 525 526 Output Parameter: 527 . A - the matrix 528 529 Level: intermediate 530 531 Notes: This matrix object does not support most matrix operations, include 532 MatSetValues(). 533 You must NOT free the ii, values and jj arrays yourself. PETSc will free them 534 when the matrix is destroyed; you must allocate them with PetscMalloc(). If you 535 call from Fortran you need not create the arrays with PetscMalloc(). 536 Should not include the matrix diagonals. 537 538 If you already have a matrix, you can create its adjacency matrix by a call 539 to MatConvert, specifying a type of MATMPIADJ. 540 541 Possible values for MatSetOption() - MAT_STRUCTURALLY_SYMMETRIC 542 543 .seealso: MatCreate(), MatConvert(), MatGetOrdering() 544 @*/ 545 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateMPIAdj(MPI_Comm comm,PetscInt m,PetscInt N,PetscInt *i,PetscInt *j,PetscInt *values,Mat *A) 546 { 547 PetscErrorCode ierr; 548 549 PetscFunctionBegin; 550 ierr = MatCreate(comm,A);CHKERRQ(ierr); 551 ierr = MatSetSizes(*A,m,PETSC_DETERMINE,PETSC_DETERMINE,N);CHKERRQ(ierr); 552 ierr = MatSetType(*A,MATMPIADJ);CHKERRQ(ierr); 553 ierr = MatMPIAdjSetPreallocation(*A,i,j,values);CHKERRQ(ierr); 554 PetscFunctionReturn(0); 555 } 556 557 558 559 560 561 562 563