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)!=NULL) { 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 with length 2n+2. First n+1 entries are set based on the traditional layout of L and U matrices 60 - bdiag - int array holding the number of nonzeros in each row of L matrix, excluding diagonals. 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_newdatastruct() for detailed description. 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; 74 PetscInt nnzL,nnzU; 75 76 PetscFunctionBegin; 77 bi[2*n+1] = 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 (bi[row+1] <= total && row < n){ 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[2*n-row+1]-1; 102 space[bdiag[row]] = row; 103 104 /* U part */ 105 nnzU = nnz - nnzL; 106 bi[2*n-row] = bi[2*n-row+1] - nnzU; 107 nnzU --; /* exclude diagonal */ 108 bj = space + bi[2*n-(row)]; 109 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 110 111 array += nnz; 112 row++; 113 } 114 115 a = (*head)->more_space; 116 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 117 ierr = PetscFree(*head);CHKERRQ(ierr); 118 *head = a; 119 } 120 bi[n] = bi[n-1] + nnzL; 121 if (bi[n] != bi[n+1]) SETERRQ2(1,"bi[n] %d != bi[n+1] %d",bi[n],bi[n+1]); 122 PetscFunctionReturn(0); 123 } 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "PetscFreeSpaceContiguous_LU_v2" 127 PetscErrorCode PetscFreeSpaceContiguous_LU_v2(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *bi,PetscInt *bdiag) 128 { 129 PetscFreeSpaceList a; 130 PetscErrorCode ierr; 131 PetscInt row,nnz,*bj,*array,total,bi_temp; 132 PetscInt nnzL,nnzU; 133 134 PetscFunctionBegin; 135 /* bi[2*n+1] = bi[n]; */ 136 bi_temp = bi[n]; 137 row = 0; 138 total = 0; 139 nnzL = bdiag[0]; 140 while ((*head)!=NULL) { 141 total += (*head)->local_used; 142 array = (*head)->array_head; 143 144 while (bi[row+1] <= total && row < n){ 145 /* copy array entries into bj for this row */ 146 nnz = bi[row+1] - bi[row]; 147 /* set bi[row] for new datastruct */ 148 if (row == 0 ){ 149 bi[row] = 0; 150 } else { 151 bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */ 152 } 153 154 /* L part */ 155 nnzL = bdiag[row]; 156 bj = space+bi[row]; 157 ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr); 158 159 /* diagonal entry */ 160 /* bdiag[row] = bi[2*n-row+1]-1; */ 161 bdiag[row] = bi_temp - 1; 162 space[bdiag[row]] = row; 163 164 /* U part */ 165 nnzU = nnz - nnzL; 166 /* bi[2*n-row] = bi[2*n-row+1] - nnzU; */ 167 bi_temp = bi_temp - nnzU; 168 nnzU --; /* exclude diagonal */ 169 /* bj = space + bi[2*n-(row)]; */ 170 bj = space + bi_temp; 171 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 172 173 array += nnz; 174 row++; 175 } 176 177 a = (*head)->more_space; 178 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 179 ierr = PetscFree(*head);CHKERRQ(ierr); 180 *head = a; 181 } 182 bi[n] = bi[n-1] + nnzL; 183 bdiag[n] = bdiag[n-1]-1; 184 /* if (bi[n] != bi[n+1]) SETERRQ2(1,"bi[n] %d != bi[n+1] %d",bi[n],bi[n+1]); */ 185 if(bi[n] != bdiag[n-1]) SETERRQ2(1,"bi[n] %d != bdiag[n-1] %d",bi[n],bdiag[n-1]); 186 PetscFunctionReturn(0); 187 } 188 189 /* 190 PetscFreeSpaceContiguous_Cholesky - 191 Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array 192 that enables an efficient matrix triangular solve. 193 194 Input Parameters: 195 + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization 196 . space - an allocated int array with length nnz of factored matrix. 197 . n - order of the matrix 198 . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix. 199 - udiag - int array of length n. 200 201 Output Parameter: 202 + space - column indices are copied into this int array with contiguous layout of U, with diagonal located as the last entry in each row 203 - udiag - indices of diagonal entries 204 205 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description. 206 */ 207 208 #undef __FUNCT__ 209 #define __FUNCT__ "PetscFreeSpaceContiguous_Cholesky" 210 PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag) 211 { 212 PetscFreeSpaceList a; 213 PetscErrorCode ierr; 214 PetscInt row,nnz,*uj,*array,total; 215 216 PetscFunctionBegin; 217 row = 0; 218 total = 0; 219 while ((*head)!=NULL) { 220 total += (*head)->local_used; 221 array = (*head)->array_head; 222 223 while (ui[row+1] <= total && row < n){ 224 udiag[row] = ui[row+1] - 1; /* points to the last entry of U(row,:) */ 225 nnz = ui[row+1] - ui[row] - 1; /* exclude diagonal */ 226 uj = space + ui[row]; 227 ierr = PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));CHKERRQ(ierr); 228 uj[nnz] = array[0]; /* diagonal */ 229 array += nnz + 1; 230 row++; 231 } 232 233 a = (*head)->more_space; 234 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 235 ierr = PetscFree(*head);CHKERRQ(ierr); 236 *head = a; 237 } 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "PetscFreeSpaceDestroy" 243 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head) 244 { 245 PetscFreeSpaceList a; 246 PetscErrorCode ierr; 247 248 PetscFunctionBegin; 249 while ((head)!=NULL) { 250 a = (head)->more_space; 251 ierr = PetscFree((head)->array_head);CHKERRQ(ierr); 252 ierr = PetscFree(head);CHKERRQ(ierr); 253 head = a; 254 } 255 PetscFunctionReturn(0); 256 } 257