1 #define PETSCMAT_DLL 2 3 #include "../src/mat/utils/freespace.h" 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "PetscFreeSpaceGet" 7 PetscErrorCode PetscFreeSpaceGet(PetscInt n,PetscFreeSpaceList *list) 8 { 9 PetscFreeSpaceList a; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = PetscMalloc(sizeof(struct _Space),&a);CHKERRQ(ierr); 14 ierr = PetscMalloc(n*sizeof(PetscInt),&(a->array_head));CHKERRQ(ierr); 15 a->array = a->array_head; 16 a->local_remaining = n; 17 a->local_used = 0; 18 a->total_array_size = 0; 19 a->more_space = NULL; 20 21 if (*list) { 22 (*list)->more_space = a; 23 a->total_array_size = (*list)->total_array_size; 24 } 25 26 a->total_array_size += n; 27 *list = a; 28 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "PetscFreeSpaceContiguous" 33 PetscErrorCode PetscFreeSpaceContiguous(PetscFreeSpaceList *head,PetscInt *space) 34 { 35 PetscFreeSpaceList a; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 while ((*head)) { 40 a = (*head)->more_space; 41 ierr = PetscMemcpy(space,(*head)->array_head,((*head)->local_used)*sizeof(PetscInt));CHKERRQ(ierr); 42 space += (*head)->local_used; 43 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 44 ierr = PetscFree(*head);CHKERRQ(ierr); 45 *head = a; 46 } 47 PetscFunctionReturn(0); 48 } 49 50 /* 51 PetscFreeSpaceContiguous_LU - 52 Copy a linket list obtained from matrix symbolic ILU or LU factorization into a contiguous array 53 that enables an efficient matrix triangular solve. 54 55 Input Parameters: 56 + head - linked list of column indices obtained from matrix symbolic ILU or LU factorization 57 . space - an allocated int array with length nnz of factored matrix. 58 . n - order of the matrix 59 . bi - row pointer of factored matrix L with length n+1. 60 - bdiag - int array of length n+1. bdiag[i] points to diagonal of U(i,:), and bdiag[n] points to entry of U(n-1,0)-1. 61 62 Output Parameter: 63 . space - column indices are copied into this int array with contiguous layout of L and U 64 65 See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed data structure of L and U 66 */ 67 #undef __FUNCT__ 68 #define __FUNCT__ "PetscFreeSpaceContiguous_LU" 69 PetscErrorCode PetscFreeSpaceContiguous_LU(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag) 70 { 71 PetscFreeSpaceList a; 72 PetscErrorCode ierr; 73 PetscInt row,nnz,*bj,*array,total,bi_temp; 74 PetscInt nnzL,nnzU; 75 76 PetscFunctionBegin; 77 bi_temp = bi[n]; 78 row = 0; 79 total = 0; 80 nnzL = bdiag[0]; 81 while ((*head)!=NULL) { 82 total += (*head)->local_used; 83 array = (*head)->array_head; 84 85 while (row < n && bi[row+1] <= total) { 86 /* copy array entries into bj for this row */ 87 nnz = bi[row+1] - bi[row]; 88 /* set bi[row] for new datastruct */ 89 if (row == 0 ){ 90 bi[row] = 0; 91 } else { 92 bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */ 93 } 94 95 /* L part */ 96 nnzL = bdiag[row]; 97 bj = space+bi[row]; 98 ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr); 99 100 /* diagonal entry */ 101 bdiag[row] = bi_temp - 1; 102 space[bdiag[row]] = row; 103 104 /* U part */ 105 nnzU = nnz - nnzL; 106 bi_temp = bi_temp - nnzU; 107 nnzU --; /* exclude diagonal */ 108 bj = space + bi_temp; 109 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 110 array += nnz; 111 row++; 112 } 113 114 a = (*head)->more_space; 115 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 116 ierr = PetscFree(*head);CHKERRQ(ierr); 117 *head = a; 118 } 119 bi[n] = bi[n-1] + nnzL; 120 bdiag[n] = bdiag[n-1]-1; 121 PetscFunctionReturn(0); 122 } 123 124 /* 125 PetscFreeSpaceContiguous_Cholesky - 126 Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array 127 that enables an efficient matrix triangular solve. 128 129 Input Parameters: 130 + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization 131 . space - an allocated int array with length nnz of factored matrix. 132 . n - order of the matrix 133 . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix. 134 - udiag - int array of length n. 135 136 Output Parameter: 137 + space - column indices are copied into this int array with contiguous layout of U, with diagonal located as the last entry in each row 138 - udiag - indices of diagonal entries 139 140 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description. 141 */ 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "PetscFreeSpaceContiguous_Cholesky" 145 PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag) 146 { 147 PetscFreeSpaceList a; 148 PetscErrorCode ierr; 149 PetscInt row,nnz,*uj,*array,total; 150 151 PetscFunctionBegin; 152 row = 0; 153 total = 0; 154 while ((*head)!=NULL) { 155 total += (*head)->local_used; 156 array = (*head)->array_head; 157 158 while (ui[row+1] <= total && row < n){ 159 udiag[row] = ui[row+1] - 1; /* points to the last entry of U(row,:) */ 160 nnz = ui[row+1] - ui[row] - 1; /* exclude diagonal */ 161 uj = space + ui[row]; 162 ierr = PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));CHKERRQ(ierr); 163 uj[nnz] = array[0]; /* diagonal */ 164 array += nnz + 1; 165 row++; 166 } 167 168 a = (*head)->more_space; 169 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 170 ierr = PetscFree(*head);CHKERRQ(ierr); 171 *head = a; 172 } 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "PetscFreeSpaceDestroy" 178 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head) 179 { 180 PetscFreeSpaceList a; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 while ((head)!=NULL) { 185 a = (head)->more_space; 186 ierr = PetscFree((head)->array_head);CHKERRQ(ierr); 187 ierr = PetscFree(head);CHKERRQ(ierr); 188 head = a; 189 } 190 PetscFunctionReturn(0); 191 } 192