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