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) { 86 if (bi[row+1] > total) break; 87 /* copy array entries into bj for this row */ 88 nnz = bi[row+1] - bi[row]; 89 /* set bi[row] for new datastruct */ 90 if (row == 0 ){ 91 bi[row] = 0; 92 } else { 93 bi[row] = bi[row-1] + nnzL; /* nnzL of previous row */ 94 } 95 96 /* L part */ 97 nnzL = bdiag[row]; 98 bj = space+bi[row]; 99 ierr = PetscMemcpy(bj,array,nnzL*sizeof(PetscInt));CHKERRQ(ierr); 100 101 /* diagonal entry */ 102 bdiag[row] = bi_temp - 1; 103 space[bdiag[row]] = row; 104 105 /* U part */ 106 nnzU = nnz - nnzL; 107 bi_temp = bi_temp - nnzU; 108 nnzU --; /* exclude diagonal */ 109 bj = space + bi_temp; 110 ierr = PetscMemcpy(bj,array+nnzL+1,nnzU*sizeof(PetscInt));CHKERRQ(ierr); 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 bdiag[n] = bdiag[n-1]-1; 122 PetscFunctionReturn(0); 123 } 124 125 /* 126 PetscFreeSpaceContiguous_Cholesky - 127 Copy a linket list obtained from matrix symbolic ICC or Cholesky factorization into a contiguous array 128 that enables an efficient matrix triangular solve. 129 130 Input Parameters: 131 + head - linked list of column indices obtained from matrix symbolic ICC or Cholesky factorization 132 . space - an allocated int array with length nnz of factored matrix. 133 . n - order of the matrix 134 . ui - row pointer of factored matrix with length n+1. All entries are set based on the traditional layout U matrix. 135 - udiag - int array of length n. 136 137 Output Parameter: 138 + space - column indices are copied into this int array with contiguous layout of U, with diagonal located as the last entry in each row 139 - udiag - indices of diagonal entries 140 141 See MatICCFactorSymbolic_SeqAIJ_newdatastruct() for detailed description. 142 */ 143 144 #undef __FUNCT__ 145 #define __FUNCT__ "PetscFreeSpaceContiguous_Cholesky" 146 PetscErrorCode PetscFreeSpaceContiguous_Cholesky(PetscFreeSpaceList *head,PetscInt *space,PetscInt n,PetscInt *ui,PetscInt *udiag) 147 { 148 PetscFreeSpaceList a; 149 PetscErrorCode ierr; 150 PetscInt row,nnz,*uj,*array,total; 151 152 PetscFunctionBegin; 153 row = 0; 154 total = 0; 155 while (*head) { 156 total += (*head)->local_used; 157 array = (*head)->array_head; 158 159 while (row < n){ 160 if (ui[row+1] > total) break; 161 udiag[row] = ui[row+1] - 1; /* points to the last entry of U(row,:) */ 162 nnz = ui[row+1] - ui[row] - 1; /* exclude diagonal */ 163 uj = space + ui[row]; 164 ierr = PetscMemcpy(uj,array+1,nnz*sizeof(PetscInt));CHKERRQ(ierr); 165 uj[nnz] = array[0]; /* diagonal */ 166 array += nnz + 1; 167 row++; 168 } 169 170 a = (*head)->more_space; 171 ierr = PetscFree((*head)->array_head);CHKERRQ(ierr); 172 ierr = PetscFree(*head);CHKERRQ(ierr); 173 *head = a; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "PetscFreeSpaceDestroy" 180 PetscErrorCode PetscFreeSpaceDestroy(PetscFreeSpaceList head) 181 { 182 PetscFreeSpaceList a; 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 while ((head)!=NULL) { 187 a = (head)->more_space; 188 ierr = PetscFree((head)->array_head);CHKERRQ(ierr); 189 ierr = PetscFree(head);CHKERRQ(ierr); 190 head = a; 191 } 192 PetscFunctionReturn(0); 193 } 194