1 /* # ***** sort to /src/direct
2 # ***** */
3 #include "mulGlobal.h"
4
Q2PDiag(chgs,numchgs,is_dummy,calc)5 double **Q2PDiag(chgs, numchgs, is_dummy, calc)
6 charge **chgs;
7 int numchgs, *is_dummy;
8 int calc;
9 {
10 double **mat;
11 int i, j;
12 double calcp();
13
14 /* Allocate storage for the potential coefficients. */
15 CALLOC(mat, numchgs, double*, ON, AQ2PD);
16 for(i=0; i < numchgs; i++) CALLOC(mat[i], numchgs, double, ON, AQ2PD);
17
18 if(calc) {
19 /* Compute the potential coeffs. */
20 /* - exclude dummy panels when they would need to carry charge
21 - exclude dielec i/f panels when they would lead to evals at their
22 centers (only if using two-point flux-den-diff evaluations) */
23 for(i=0; i < numchgs; i++) {
24 #if NUMDPT == 2
25 if(chgs[i]->dummy); /* don't check surface of a dummy */
26 else if(chgs[i]->surf->type == DIELEC || chgs[i]->surf->type == BOTH)
27 continue;
28 #endif
29 for(j=0; j < numchgs; j++) { /* need to have charge on them */
30 #if SKIPQD == ON
31 if(chgs[j]->pos_dummy == chgs[i] || chgs[j]->neg_dummy == chgs[i])
32 continue;
33 #endif
34 if(!is_dummy[j]) mat[i][j] = calcp(chgs[j], chgs[i], NULL);
35 }
36 }
37 }
38
39 #if DSQ2PD == ON
40 dispQ2PDiag(mat, chgs, numchgs, is_dummy);
41 #endif
42
43 return(mat);
44 }
45
Q2P(qchgs,numqchgs,is_dummy,pchgs,numpchgs,calc)46 double **Q2P(qchgs, numqchgs, is_dummy, pchgs, numpchgs, calc)
47 charge **qchgs, **pchgs;
48 int numqchgs, numpchgs;
49 int *is_dummy;
50 int calc;
51 {
52 double **mat;
53 int i, j;
54 double calcp();
55
56 /* Allocate storage for the potential coefficients. P rows by Q cols. */
57 CALLOC(mat, numpchgs, double*, ON, AQ2P);
58 for(i=0; i < numpchgs; i++) {
59 CALLOC(mat[i], numqchgs, double, ON, AQ2P);
60 if(calc) {
61 /* exclude:
62 - dummy panels in the charge list
63 - dielectric i/f panels in the eval list (if doing 2-point E's)*/
64 #if NUMDPT == 2
65 if(pchgs[i]->dummy); /* don't check the surface of a dummy */
66 else if(pchgs[i]->surf->type == DIELEC || pchgs[i]->surf->type == BOTH)
67 continue;
68 #endif
69 for(j=0; j < numqchgs; j++) { /* only dummy panels in the charge list */
70 if(!is_dummy[j]) /* (not the eval list) are excluded */
71 mat[i][j] = calcp(qchgs[j], pchgs[i], NULL);
72 /* old: pchgs[i],qchgs[j] */
73 }
74 }
75 }
76
77 #if DISQ2P == ON
78 dispQ2P(mat, qchgs, numqchgs, is_dummy, pchgs, numpchgs);
79 #endif
80
81 return(mat);
82 }
83
84 /*
85 used only in conjunction with DMPMAT == ON and DIRSOL == ON
86 to make 1st directlist mat = full P mat
87 */
Q2Pfull(directlist,numchgs)88 double **Q2Pfull(directlist, numchgs)
89 int numchgs;
90 cube *directlist;
91 {
92 int i, j, fromp, fromq, top, toq;
93 double **mat, calcp();
94 cube *pq, *pp;
95 charge **pchgs, **qchgs, *eval;
96
97 /* allocate room for full P matrix */
98 MALLOC(mat, numchgs, double*, ON, AQ2P);
99 for(i = 0; i < numchgs; i++) MALLOC(mat[i], numchgs, double, ON, AQ2P);
100
101 /* load the matrix in the style of Q2P() - no attempt to exploit symmetry */
102 /* calculate interaction between every direct list entry and entire dlist */
103 for(pp = directlist; pp != NULL; pp = pp->dnext) {
104 pchgs = pp->chgs;
105 fromp = pchgs[0]->index - 1; /* row index range */
106 top = fromp + pp->upnumeles[0];
107 for(pq = directlist; pq != NULL; pq = pq->dnext) {
108 qchgs = pq->chgs;
109 fromq = qchgs[0]->index - 1; /* column index range */
110 toq = fromq + pq->upnumeles[0];
111
112 for(i = fromp; i < top; i++) {
113 for(j = fromq; j < toq; j++) {
114 eval = qchgs[j-fromq];
115 mat[i][j] = calcp(pchgs[i-fromp],eval, NULL);
116 }
117 }
118
119 }
120 }
121 return(mat);
122 }
123
124 /*
125 - returned matrix has L below the diagonal, U above (GVL1 pg 58)
126 - if allocate == TRUE ends up storing P and LU (could be a lot)
127 */
ludecomp(matin,size,allocate)128 double **ludecomp(matin, size, allocate)
129 double **matin;
130 int size;
131 int allocate;
132 {
133 extern int fulldirops;
134 double factor, **mat;
135 int i, j, k;
136
137 if(allocate == TRUE) {
138 /* allocate for LU matrix and copy A */
139 MALLOC(mat, size, double*, ON, AMSC);
140 for(i = 0; i < size; i++) {
141 MALLOC(mat[i], size, double, ON, AMSC);
142 for(j = 0; j < size; j++) mat[i][j] = matin[i][j];
143 }
144 }
145 else mat = matin;
146
147 for(k = 0; k < size-1; k++) { /* loop on rows */
148 if(mat[k][k] == 0.0) {
149 fprintf(stderr, "ludecomp: zero piovt\n");
150 exit(1);
151 }
152 for(i = k+1; i < size; i++) { /* loop on remaining rows */
153 factor = (mat[i][k] /= mat[k][k]);
154 fulldirops++;
155 for(j = k+1; j < size; j++) { /* loop on remaining columns */
156 mat[i][j] -= (factor*mat[k][j]);
157 fulldirops++;
158 }
159 }
160 }
161 return(mat);
162 }
163
164 /*
165 For direct solution of Pq = psi, used if DIRSOL == ON or if preconditioning.
166 */
solve(mat,x,b,size)167 void solve(mat, x, b, size)
168 double **mat, *x, *b;
169 int size;
170 {
171 extern int fulldirops;
172 int i, j;
173
174 /* copy rhs */
175 if(x != b) for(i = 0; i < size; i++) x[i] = b[i];
176
177 /* forward elimination */
178 for(i = 0; i < size; i++) { /* loop on pivot row */
179 for(j = i+1; j < size; j++) { /* loop on elimnation row */
180 x[j] -= mat[j][i]*x[i];
181 fulldirops++;
182 }
183 }
184
185 /* back substitution */
186 for(i--; i > -1; i--) { /* loop on rows */
187 for(j = i+1; j < size; j++) { /* loop on columns */
188 x[i] -= mat[i][j]*x[j];
189 fulldirops++;
190 }
191 x[i] /= mat[i][i];
192 fulldirops++;
193 }
194 }
195
196 /*
197 In-place inverts a matrix using guass-jordan.
198 - is_dummy[i] = 0 => ignore row/col i
199 */
invert(mat,size,reorder)200 void invert(mat, size, reorder)
201 double **mat;
202 int size, *reorder;
203 {
204 int i, j, k, best;
205 double normal, multiplier, bestval, nextbest;
206 /*
207 matlabDump(mat,size,"p");
208 */
209 for(i=0; i < size; i++) {
210 best = i;
211 bestval = ABS(mat[i][i]);
212 for(j = i+1; j < size; j++) {
213 nextbest = ABS(mat[i][j]);
214 if(nextbest > bestval) {
215 best = j;
216 bestval = nextbest;
217 }
218 }
219
220 /* If reordering, find the best pivot. */
221 if(reorder != NULL) {
222 reorder[i] = best;
223 if(best != i) {
224 for(k=0; k < size; k++) {
225 bestval = mat[k][best];
226 mat[k][best] = mat[k][i];
227 mat[k][i] = bestval;
228 }
229 }
230 }
231
232 /* First i^{th} column of A. */
233 normal = 1.0 / mat[i][i];
234 for(j=0; j < size; j++) {
235 mat[j][i] *= normal;
236 }
237 mat[i][i] = normal;
238
239 /* Fix the backward columns. */
240 for(j=0; j < size; j++) {
241 if(j != i) {
242 multiplier = -mat[i][j];
243 for(k=0; k < size; k++) {
244 if(k != i) mat[k][j] += mat[k][i] * multiplier;
245 else mat[k][j] = mat[k][i] * multiplier;
246 }
247 }
248 }
249 }
250
251 /* Unravel the reordering, starting with the last column. */
252 if(reorder != NULL) {
253 for(i=size-2; i >= 0; i--) {
254 if(reorder[i] != i) {
255 for(k=0; k < size; k++) {
256 bestval = mat[k][i];
257 mat[k][reorder[i]] = mat[k][i];
258 mat[k][i] = bestval;
259 }
260 }
261 }
262 }
263 /*
264 matlabDump(mat,size,"c");
265 */
266
267 }
268
269 /*
270 Used in conjuction with invert() to remove dummy row/columns
271 comp_rows = TRUE => remove rows corresponding to ones in is_dummy[]
272 comp_rows = FALSE => remove cols corresponding to ones in is_dummy[]
273 comp_rows = BOTH => remove both rows and columns
274 returns number of rows/cols in compressed matrix
275 */
compressMat(mat,size,is_dummy,comp_rows)276 int compressMat(mat, size, is_dummy, comp_rows)
277 double **mat;
278 int size, *is_dummy, comp_rows;
279 {
280 static int *cur_order;
281 static int cur_order_array_size = 0;
282 int cur_order_size, i, j, k;
283
284 if(cur_order_array_size < size) {
285 CALLOC(cur_order, size, int, ON, AMSC);
286 }
287
288 /* figure the new order vector (cur_order[i] = index of ith row/col) */
289 for(i = cur_order_size = 0; i < size; i++) {
290 if(!is_dummy[i]) cur_order[cur_order_size++] = i;
291 }
292
293 if(comp_rows == TRUE || comp_rows == BOTH) {
294 /* compress by removing rows from the matrix */
295 for(i = 0; i < cur_order_size; i++) {
296 if((k = cur_order[i]) == i) continue; /* if not reordered */
297 for(j = 0; j < size; j++) { /* copy the row to its new location */
298 mat[i][j] = mat[k][j];
299 }
300 }
301 }
302 if(comp_rows == FALSE || comp_rows == BOTH) {
303 /* compress by removing columns from the matrix */
304 for(j = 0; j < cur_order_size; j++) {
305 if((k = cur_order[j]) == j) continue; /* if not reordered */
306 for(i = 0; i < size; i++) { /* copy the col to its new location */
307 mat[i][j] = mat[i][k];
308 }
309 }
310 }
311 return(cur_order_size);
312 }
313
314 /*
315 Used in conjuction with invert() to add dummy row/columns
316 exp_rows = TRUE => add rows corresponding to ones in is_dummy[]
317 exp_rows = FALSE => add cols corresponding to ones in is_dummy[]
318 exp_rows = BOTH => add rows and columns
319 */
expandMat(mat,size,comp_size,is_dummy,exp_rows)320 void expandMat(mat, size, comp_size, is_dummy, exp_rows)
321 double **mat;
322 int size, *is_dummy, exp_rows, comp_size;
323 {
324 int i, j, k, next_rc;
325
326 if(exp_rows == TRUE || exp_rows == BOTH) {
327 next_rc = comp_size - 1;
328 /* add rows to the matrix starting from the bottom */
329 for(i = size - 1; i >= 0; i--) {
330 if(is_dummy[i]) { /* zero out dummy row */
331 for(j = 0; j < size; j++) mat[i][j] = 0.0;
332 }
333 else { /* copy the row from its compressed location */
334 for(j = 0; j < size; j++) mat[i][j] = mat[next_rc][j];
335 next_rc--;
336 }
337 }
338 }
339 if(exp_rows == FALSE || exp_rows == BOTH) {
340 next_rc = comp_size - 1;
341 /* add columns to the matrix starting from the right */
342 for(j = size - 1; j >= 0; j--) {
343 if(is_dummy[j]) { /* zero out dummy column */
344 for(i = 0; i < size; i++) mat[i][j] = 0.0;
345 }
346 else { /* copy the col from its compressed location */
347 for(i = 0; i < size; i++) mat[i][j] = mat[i][next_rc];
348 next_rc--;
349 }
350 }
351 }
352 }
353
354 /*
355 Checks to see if the matrix has the M-matrix sign pattern and if
356 it is diagonally dominant.
357 */
matcheck(mat,rows,size)358 matcheck(mat, rows, size)
359 double **mat;
360 int rows, size;
361 {
362 double rowsum;
363 int i, j;
364
365 for(i = rows - 1; i >= 0; i--) {
366 for(rowsum = 0.0, j = size - 1; j >= 0; j--) {
367 if((i != j) && (mat[i][j] > 0.0)) {
368 printf("violation mat[%d][%d] =%g\n", i, j, mat[i][j]);
369 }
370 if(i != j) rowsum += ABS(mat[i][j]);
371 }
372 printf("row %d diag=%g rowsum=%g\n", i, mat[i][i], rowsum);
373 if(rowsum > mat[i][i]) {
374 for(j = size - 1; j >= 0; j--) {
375 printf("col%d = %g ", j, mat[i][j]);
376 }
377 printf("\n");
378 }
379 }
380 }
381
382
matlabDump(mat,size,name)383 matlabDump(mat, size, name)
384 double **mat;
385 int size;
386 char *name;
387 {
388 FILE *foo;
389 int i,j;
390 char fname[100];
391
392 sprintf(fname, "%s.m", name);
393 /* SRW -- this is ascii data */
394 foo = fopen(fname, "w");
395 fprintf(foo, "%s = [\n", name);
396 for(i=0; i < size; i++) {
397 for(j=0; j < size; j++) {
398 fprintf(foo, "%.10e ", mat[i][j]);
399 }
400 fprintf(foo, "\n");
401 }
402 fprintf(foo, "]\n");
403 }
404
405
406
407
408
409
410
411
412