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