1 /* This is the main part of the code */
2
3 #include "induct.h"
4 #include <string.h>
5
6 /* these are missing in some math.h files */
7 extern double asinh();
8 extern double atanh();
9
10 #define MAXPRINT 1
11 #define MAXCHARS 400
12 #define TIMESIZE 10
13
14 static int notblankline();
15
16 FILE *fp, *fp2, *fp3, *fptemp, *fb, *fROM;
17 int num_exact_mutual;
18 int num_fourfil;
19 int num_mutualfil;
20 int num_found;
21 int num_perp;
22 int forced = 0; /* for debugging inside exact_mutual() */
23
24 char outfname[200];
25 char outfname2[200];
26
27 /* savemat_mod machine type */
28 #ifdef DEC
29 int machine = 0000;
30 #else
31 int machine = 1000;
32 #endif
33
34 #if SUPERCON == ON
35 void fillZ_diag(SYS *indsys, double omega);
36 void set_rvals(SYS *indsys, double omega);
37 #endif
38
39
40 int
main(argc,argv)41 main(argc, argv)
42 int argc;
43 char *argv[];
44 {
45
46 double width, height, length, freq, freqlast;
47 int Linc, Winc, Hinc, filnum;
48 double ratio;
49 GROUNDPLANE *plane; /* CMS 7/7/92 */
50 FILAMENT *tmpf;
51 SEGMENT *tmps;
52 SEGMENT *seg;
53 NODES *node;
54 int i,j,k,m, last, err;
55 char fname[80], tempstr[10];
56 double r_height, r_width;
57 CX dumb;
58 CX *vect, **pvect; /* space needed by gmres */
59 double tol = 1e-8;
60 double ftimes[TIMESIZE];
61 int maxiters, user_maxiters;
62 CX *b, *x0;
63 SYS *indsys; /* holds all the big variables (inductance system) */
64 CX **MtZM;
65 int num_mesh, num_extern, num_fils,num_nodes, num_real_nodes, num_segs;
66 int num_sub_extern;
67 double totaltime;
68 int dont_form_Z; /* if fmin=0, don't from the L matrix */
69 char *MRMt; /* may be needed if ROM is requested */
70 double **MRMtinvMLMt, **B, **C;
71 double **romB, **romC;
72 int actual_order;
73
74 int num_planes, nonp, planemeshes, tree_meshes; /* CMS 6/7/92 */
75 void fillgrids(); /* function in addgroundplane.c */
76
77 int *meshsect;
78 double fmin, fmax, logofstep;
79 double *R;
80 EXTERNAL *ext;
81 ind_opts *opts;
82 extern long memcount;
83
84 charge *chglist, *chgend, chgdummy;
85 ssystem *sys, *SetupMulti();
86
87 memcount = 0;
88
89 for(i = 0; i < TIMESIZE; i++)
90 ftimes[i] = 0;
91
92 starttimer;
93
94 indsys = (SYS *)MattAlloc(1, sizeof(SYS));
95 indsys->externals = NULL;
96 indsys->num_extern = 0;
97 indsys->segment = NULL;
98 indsys->nodes = NULL;
99 indsys->planes = NULL; /* CMS 7/2/92 */
100 indsys->logofstep = 0;
101 indsys->opts = Parse_Command_Line(argc, argv);
102
103 /*indsys->r_height = indsys->r_width = indsys->opts->filratio;*/
104
105 /* set the type equal to that specified on command line for now */
106 indsys->precond_type = indsys->opts->precond;
107
108 /* degugging counters for calls to functions used in mutual() */
109 num_exact_mutual = 0;
110 num_fourfil = 0;
111 num_mutualfil = 0;
112 num_found = 0;
113 num_perp = 0;
114
115 opts = indsys->opts;
116
117 tol = opts->tol;
118 user_maxiters = opts->maxiters;
119
120 if (opts->fname == NULL) {
121 printf("No input file given. Reading from stdin...\n");
122 fp = stdin;
123 }
124 else if (strcmp(opts->fname,"-") == 0) {
125 printf("Reading from stdin...\n");
126 fp = stdin;
127 }
128 else {
129 /* SRW -- read ascii file */
130 fp = fopen(opts->fname, "r");
131 if (fp == NULL) {
132 printf("Couldn't open %s\n", opts->fname);
133 exit(1);
134 }
135 printf("Reading from file: %s\n",opts->fname);
136 }
137
138 /* read in geometry */
139 err = readGeom(fp, indsys);
140 if (err != 0)
141 return err;
142
143 fclose(fp);
144
145 /* fprintf(stdout,"Number of nodes before breaking up: %d\n",
146 indsys->num_nodes); */
147
148 /* regurgitate input file to stdout */
149 if (opts->regurgitate == TRUE) {
150 regurgitate(indsys);
151 }
152
153 if ((opts->makeFastCapFile & HIERARCHY)
154 && !(opts->makeFastCapFile & (SIMPLE | REFINED)))
155 /* the hierarchy has been dumped already, and nothing more to do. so quit*/
156 return 0;
157
158 /* make a file suitable for keith's postscript maker */
159
160 if (opts->makeFastCapFile & SIMPLE) {
161 fprintf(stdout, "Making simple zbuf file...");
162 fflush(stdout);
163 concat4(outfname,"zbuffile",opts->suffix,"");
164 concat4(outfname2,outfname,"_","shadings");
165 writefastcap(outfname, outfname2, indsys);
166 fprintf(stdout, "Done\n");
167 if (! (opts->makeFastCapFile & REFINED) )
168 /* we are done, after making fastcap files for visualization */
169 return 0;
170 }
171
172 /* initialize Gaussian quadrature arrays for mutual(). */
173 /* See dist_betw_fils.c */
174 fill_Gquad();
175
176 /* initialize Freeable allocation stuff for mutual terms lookup table */
177 init_table();
178
179 /* break each segment into filaments */
180 num_fils = 0;
181 chgend = &chgdummy;
182
183 for(seg = indsys->segment; seg != NULL; seg = seg->next)
184 chgend = assignFil(seg, &num_fils, chgend);
185
186 indsys->num_fils = num_fils;
187
188 chgend->next = NULL;
189 chglist = chgdummy.next;
190
191 stoptimer;
192 ftimes[0] = dtime;
193
194 starttimer;
195
196 #if SUPERCON == ON
197 set_rvals(indsys, 2*PI*indsys->fmin);
198 #endif
199
200 /* set up multipole stuff */
201 sys = SetupMulti(chglist, indsys);
202
203 #if 1 == 0
204 if (indsys->opts->debug == ON && 1 == 0) /* for debugging eval pass */
205 dump_evalcnts(sys);
206 #endif
207
208 stoptimer;
209 ftimes[5] = dtime;
210
211 if (indsys->opts->debug == ON)
212 printf("Time for Multipole Setup: %lg\n",dtime);
213
214 starttimer;
215
216 printf("Scanning graph to find fundamental circuits...\n");
217 /* find all the necessary meshes */
218 make_trees(indsys);
219
220 /* clear node->examined and seg->is_deleted */
221 clear_marks(indsys);
222
223 /* find meshes in groundplane due to holes and put them at the front
224 of the list of trees so that they are created by fillM() first so
225 that they are marked before regular ground plane meshes */
226 find_hole_meshes(indsys);
227
228 /* determine if a subset of columns is to be computed and assign
229 col_Yindex. This will happen if -x option is used. */
230 num_sub_extern = indsys->num_sub_extern
231 = pick_subset(opts->portlist, indsys);
232
233 stoptimer;
234 ftimes[6] = dtime;
235
236 /* Write to Zc.mat only if we aren't only running for visualization */
237 if ( !(opts->makeFastCapFile & (SIMPLE | REFINED))
238 && !(opts->orderROM > 0 && opts->onlyROM) ) {
239 concat4(outfname,"Zc",opts->suffix,".mat"); /* put filnames together */
240 /* SRW -- this is ascii data */
241 fp3 = fopen(outfname, "w");
242 if (fp3 == NULL) {
243 printf("couldn't open file %s\n",outfname);
244 exit(1);
245 }
246
247 for(ext = indsys->externals; ext != NULL; ext=ext->next) {
248 /* printf("Row %d : %s to %s\n",ext->Yindex,ext->source->node[0]->name,
249 ext->source->node[1]->name); */
250
251 if (ext->portname[0] == '\0')
252 fprintf(fp3,"Row %d: %s to %s\n",ext->Yindex+1,ext->name1,ext->name2);
253 else
254 fprintf(fp3,"Row %d: %s to %s, port name: %s\n",
255 ext->Yindex+1,ext->name1,ext->name2,ext->portname);
256 }
257
258 if (num_sub_extern != indsys->num_extern)
259 for(ext = indsys->externals; ext != NULL; ext=ext->next) {
260 if (ext->col_Yindex != -1)
261 fprintf(fp3,"Col %d: port name: %s\n",
262 ext->col_Yindex+1,ext->portname);
263 }
264 }
265
266 /* count nodes */
267 indsys->num_real_nodes = 0;
268 indsys->num_nodes = 0;
269 last = -1;
270 for(node = indsys->nodes, i = 0; node != NULL; node = node->next, i++) {
271 indsys->num_nodes++;
272 if (getrealnode(node) == node) {
273 indsys->num_real_nodes++;
274 if (last + 1 != i) {
275 /* printf("Non equivalent nodes must be listed first??,
276 no I take that back.\n"); */
277 /* exit(1); */
278 }
279 last = i;
280 }
281 }
282
283 /* CMS 7/3/92 ------------------------------------------------------------*/
284 planemeshes = 0;
285 nonp = 0;
286
287 for(plane = indsys->planes; plane != NULL; plane = plane->next){
288 planemeshes = planemeshes + plane->numesh;
289 if(plane->external == 0)
290 nonp++;
291 }
292 /*------------------------------------------------------------------------*/
293
294 /* moved lower for shading
295 if (opts->makeFastCapFile & REFINED) {
296 fprintf(stdout, "Making refined zbuf file...");
297 fflush(stdout);
298
299 concat4(outfname,"zbuffile2",opts->suffix,"");
300
301
302 concat4(outfname2,outfname,"_","shadings");
303 writefastcap(outfname, outfname2, indsys);
304 fprintf(stdout, "Done\n");
305 }
306 */
307
308 num_fils = indsys->num_fils; /* setup multi may alter this */
309 num_segs = indsys->num_segs;
310 num_nodes = indsys->num_nodes;
311 num_real_nodes = indsys->num_real_nodes;
312 num_planes = indsys->num_planes; /* CMS 7/2/92 */
313
314 /* figure out number of meshes */
315 tree_meshes = indsys->tree_meshes = count_tree_meshes(indsys->trees);
316
317 indsys->extra_meshes = tree_meshes;
318
319 #if 1==0
320 unimplemented junk
321 indsys->extra_meshes = estimate_extra_meshes(indsys->trees, FILS_PER_MESH);
322 #endif
323
324 num_mesh = num_fils - num_segs + indsys->extra_meshes + planemeshes;
325 /* CMS 7/2/92 */
326 indsys->num_mesh = num_mesh;
327 num_extern = count_externals(indsys->externals);
328 if (num_extern != indsys->num_extern) {
329 fprintf(stderr, "main: discrepancy in num_extern and actual number\n");
330 exit(1);
331 }
332 fmin = indsys->fmin;
333 fmax = indsys->fmax;
334 logofstep = indsys->logofstep;
335
336 /* dont_form_Z = fmin == 0 && opts->mat_vect_prod == DIRECT
337 && opts->soln_technique == ITERATIVE; */
338
339 dont_form_Z = fmin == 0 && opts->soln_technique == LUDECOMP;
340 indsys->dont_form_Z = dont_form_Z;
341
342
343 printf("Number of Groundplanes : %d \n",num_planes); /* CMS 7/2/92 */
344 printf("Number of filaments: %10d\n", num_fils);
345 printf("Number of segments: %10d\n", num_segs);
346 printf("Number of nodes: %10d\n", num_nodes);
347 printf("Number of meshes: %10d\n", num_mesh);
348 printf(" ----from tree: %d \n",tree_meshes);
349 printf(" ----from planes: (before holes) %d \n",planemeshes);
350 /* CMS 7/6/92 */
351 printf("Number of conductors:%10d (rows of matrix in Zc.mat) \n",
352 num_extern);
353 printf("Number of columns: %10d (columns of matrix in Zc.mat) \n",
354 num_sub_extern);
355
356 printf("Number of real nodes:%10d\n", num_real_nodes); /* non equivalent */
357
358 /* A = MatrixAlloc(num_real_nodes, num_fils, sizeof(double)); */
359 if (opts->mat_vect_prod == DIRECT || opts->soln_technique == LUDECOMP) {
360 /* indsys->M = MatrixAlloc(num_fils, num_mesh, sizeof(double)); */
361 if (!dont_form_Z)
362 indsys->Z = MatrixAlloc(num_fils, num_fils, sizeof(double));
363
364 /* let's save space and allocate this later
365 /* indsys->MtZM = (CX **)MatrixAlloc(num_mesh, num_mesh, sizeof(CX)); */
366 }
367 indsys->FinalY = (CX **)MatrixAlloc(num_extern, num_sub_extern, sizeof(CX));
368 b = (CX *)MattAlloc(num_mesh,sizeof(CX));
369 x0 = (CX *)MattAlloc(num_mesh,sizeof(CX));
370 indsys->R = (double *)MattAlloc(num_fils,sizeof(double));
371 indsys->meshsect = (int *)MattAlloc((num_extern+nonp+1),sizeof(int));
372 vect = (CX *)MattAlloc(num_mesh,sizeof(CX));
373 pvect = (CX **)MattAlloc(num_mesh,sizeof(CX *));
374 indsys->Mlist = (MELEMENT **)MattAlloc(num_mesh,sizeof(MELEMENT *));
375 indsys->Mtrans = (MELEMENT **)MattAlloc(num_fils,sizeof(MELEMENT *));
376 indsys->m_info = (Minfo *)MattAlloc(num_mesh,sizeof(Minfo));
377 indsys->Precond = (PRE_ELEMENT **)MattAlloc(num_mesh, sizeof(PRE_ELEMENT *));
378 if (1 == 1) { /* let's always dump the residual history */
379 indsys->resids = MatrixAlloc(num_sub_extern,user_maxiters, sizeof(double));
380 indsys->resid_real = MatrixAlloc(num_sub_extern, user_maxiters,
381 sizeof(double));
382 indsys->resid_imag = MatrixAlloc(num_sub_extern, user_maxiters,
383 sizeof(double));
384 indsys->niters = (double *)MattAlloc(num_sub_extern, sizeof(double));
385 }
386
387 /*
388 if (opts->mat_vect_prod == DIRECT || opts->soln_technique == LUDECOMP) {
389 for(i = 0; i < num_fils; i++) {
390 for(j = 0; j < num_mesh; j++)
391 indsys->M[i][j] = 0;
392 }
393 }
394 */
395
396 meshsect = indsys->meshsect;
397 R = indsys->R;
398
399 starttimer;
400 /*
401 printf("filling A...\n");
402 fillA(segment, A, num_segs);
403 */
404
405 printf("filling M...\n");
406 fillM(indsys); /* expects hole meshes to be marked */
407
408 if (opts->makeFastCapFile & REFINED) {
409 fprintf(stdout, "Making refined zbuf file...");
410 fflush(stdout);
411 concat4(outfname,"zbuffile2",opts->suffix,"");
412
413 /* shadings file */
414 concat4(outfname2,outfname,"_","shadings");
415
416 /* output file to later turn into .ps. Note: this affects seg->is_deleted*/
417 writefastcap(outfname, outfname2, indsys);
418 fprintf(stdout, "Done\n");
419 return 0; /* we are done after making fastcap files for visualization */;
420 }
421
422
423 if (num_mesh != indsys->num_mesh)
424 printf("Exact number of meshes after ground plane holes removed: %d\n",
425 indsys->num_mesh);
426
427 num_mesh = indsys->num_mesh; /* actual number of meshes */
428
429 /* let's save a little space and allocate MZMt now. */
430 if (!dont_form_Z
431 && (opts->mat_vect_prod == DIRECT || opts->soln_technique == LUDECOMP))
432 indsys->MtZM = (CX **)MatrixAlloc(num_mesh, num_mesh, sizeof(CX));
433
434 MtZM = indsys->MtZM;
435
436 if (opts->dumpMats & MESHES) {
437 if (opts->kind & MATLAB)
438 dump_mesh_coords(indsys);
439 else
440 dump_ascii_mesh_coords(indsys);
441 }
442
443 formMtrans(indsys); /* Form M transpose by row*/
444
445 printf("filling R and L...\n");
446 fillZ(indsys);
447
448 stoptimer;
449 ftimes[1] = dtime;
450
451 if (indsys->opts->debug == ON)
452 printf("Time to Form M and Z: %lg\n",dtime);
453
454 printf("Total Memory allocated: %ld kilobytes\n",memcount/1024);
455
456 if (indsys->opts->debug == ON)
457 printf("Memory used and freed by lookup table: %d kilobytes\n",
458 get_table_mem());
459
460 /* free memory for lookup table */
461 destroy_table();
462
463 starttimer;
464 choose_and_setup_precond(indsys);
465
466 if (opts->dumpMats) {
467 printf("saving some files disk...\n");
468 savemats(indsys);
469 }
470
471 if (logofstep == 0.0) {
472 printf("no frequency range data read!\n");
473 exit(1);
474 }
475
476 if (fmin == 0) {
477 printf("***First frequency is zero. Only the DC case will be run.***\n");
478 if (!dont_form_Z)
479 printf("Warning: First frequency is zero, but -sludecomp was not specified.\n\
480 Use this setting to save time and memory.\n");
481 }
482
483 if (indsys->opts->debug == ON) {
484 /* open Ycond.mat */
485 concat4(outfname,"Ycond",opts->suffix,".mat");
486 /* SRW -- this is binary data */
487 fp = fopen(outfname, "wb");
488 if (fp == NULL) {
489 printf("couldn't open file %s\n",outfname);
490 exit(1);
491 }
492 }
493
494 /* open b.mat */
495 if (opts->dumpMats != OFF) {
496 concat4(outfname,"b",opts->suffix,".mat");
497 /* SRW -- this is binary data */
498 if ((fb = fopen(outfname,"wb")) == NULL)
499 fprintf(stderr, "No open fb\n");
500 }
501
502 /* Model Order Reduction: create, compute and print the model if requested */
503 if (opts->orderROM > 0) {
504 double *dtemp;
505
506 /* create the matrix whose inverse we will need */
507 createMRMt(&MRMt, indsys);
508 /* now create the input and output matrices for the state */
509 B = MatrixAlloc(num_mesh, num_extern, sizeof(double));
510 C = MatrixAlloc(num_mesh, num_extern, sizeof(double));
511 for(ext = indsys->externals, j = 0;
512 ext != NULL; ext = ext->next, j++) {
513 int_list *elem;
514 /* create C first; C = N, where Vs = N Vterminal */
515 for(elem = ext->indices; elem != NULL; elem = elem->next) {
516 if (elem->index > num_mesh || ext->Yindex > num_extern) {
517 fprintf(stderr, "Indexing into matrix C out of bounds\n");
518 exit(1);
519 }
520 C[elem->index][ext->Yindex] = elem->sign;
521 B[elem->index][ext->Yindex] = elem->sign;
522 }
523 }
524
525 if (indsys->opts->mat_vect_prod != MULTIPOLE && !indsys->dont_form_Z) {
526 printf("multiplying M*(L)*transpose(M)for model order reduction\n");
527 /* this form of storage isn't the best */
528 formMLMt(indsys); /*form (M^t)*(L)*M and store in indys->MtZM*/
529
530
531 actual_order = ArnoldiROM(B, C, (double **)NULL, MRMt, num_mesh,
532 num_extern, num_extern, opts->orderROM,
533 realMatVect, indsys, sys, chglist);
534 }
535 else if (indsys->opts->mat_vect_prod == MULTIPOLE)
536 actual_order = ArnoldiROM(B, C, (double **)NULL, MRMt, num_mesh,
537 num_extern, num_extern, opts->orderROM,
538 realComputePsi, indsys, sys, chglist);
539
540 if (indsys->opts->debug == ON) {
541 #if 1==0
542 /* open orig.mat */
543 concat4(outfname,"orig",opts->suffix,".mat");
544 /* SRW -- this is binary data */
545 if ((fROM = fopen(outfname,"wb")) == NULL)
546 fprintf(stderr, "No open fROM\n");
547 /* dump what we have of the original system */
548 dumpROMbin(fROM, NULL, B, C, NULL,
549 num_mesh, num_extern, num_extern);
550 fclose(fROM);
551 #endif
552
553 /* open rom.m */
554 concat4(outfname,"rom",opts->suffix,".m");
555 /* SRW -- this is ascii data */
556 if ((fROM = fopen(outfname,"w")) == NULL)
557 fprintf(stderr, "No open fROM\n");
558
559 /* now dump the reduced order model */
560 dumpROM(fROM, indsys->Ar, indsys->Br, indsys->Cr, indsys->Dr,
561 actual_order * num_extern, num_extern, num_extern);
562 /* close and save away the file */
563 fclose(fROM);
564 }
565
566 /* generate equivalent circuit */
567 concat4(outfname,"equiv_circuitROM",opts->suffix,".spice");
568 /* SRW -- this is ascii data */
569 if ((fROM = fopen(outfname,"w")) == NULL)
570 fprintf(stderr, "No open fROM\n");
571 /* now dump the reduced order model */
572 dumpROMequiv_circuit(fROM, indsys->Ar, indsys->Br, indsys->Cr, indsys->Dr,
573 actual_order * num_extern, num_extern, num_extern,
574 indsys->title, opts->suffix, indsys);
575 /* close and save away the file */
576 fclose(fROM);
577 /* end of Model Order Reduction generation */
578 }
579
580 /* NOTE: exit here if all we wnat is the reduced order model (ROM) */
581 if (opts->orderROM > 0 && opts->onlyROM)
582 exit(0);
583
584 for(freq = fmin, m = 0;
585 (fmin != 0 && freq <= fmax*1.001) || (fmin == 0 && m == 0);
586 m++, freq = (fmin != 0 ? pow(10.0,log10(fmin) + m*logofstep) : 0.0)) {
587 printf("Frequency = %lg\n",freq);
588
589 starttimer;
590 #if SUPERCON == ON
591 if (freq != fmin)
592 fillZ_diag(indsys, 2*PI*freq);
593 #endif
594
595 if (!dont_form_Z
596 && (opts->mat_vect_prod == DIRECT || opts->soln_technique==LUDECOMP)) {
597 printf("multiplying M*(R + jL)*transpose(M)\n");
598 formMZMt(indsys); /*form transpose(M)*(R+jL)*M (no w) */
599
600 if (opts->dumpMats & MZMt) {
601 if (m == 0) {
602 if (opts->kind & MATLAB) {
603 concat4(outfname,"MZMt",opts->suffix,".mat");
604 /* SRW -- this is binary data */
605 if ( (fp2 = fopen(outfname,"wb")) == NULL) {
606 printf("Couldn't open file\n");
607 exit(1);
608 }
609 printf("Saving MZMt...\n");
610 savecmplx2(fp2,"MZMt",indsys->MtZM, indsys->num_mesh,indsys->num_mesh);
611 fclose(fp2);
612 }
613 if (opts->kind & TEXT) {
614 concat4(outfname,"MZMt",opts->suffix,".dat");
615 /* SRW -- this is ascii data */
616 if ( (fp2 = fopen(outfname,"w")) == NULL) {
617 printf("Couldn't open file\n");
618 exit(1);
619 }
620 cx_dumpMat_totextfile(fp2, indsys->MtZM,
621 indsys->num_mesh,indsys->num_mesh );
622 fclose(fp2);
623 }
624 }
625 }
626
627 printf("putting in frequency \n");
628
629 /* put in frequency */
630 for(i = 0; i < num_mesh; i++)
631 for(j = 0; j < num_mesh; j++)
632 MtZM[i][j].imag *= 2*PI*freq;
633 }
634
635 stoptimer;
636 ftimes[2] += dtime;
637
638 starttimer;
639
640 if (opts->soln_technique == ITERATIVE) {
641 if (indsys->precond_type == LOC) {
642 printf("Forming local inversion preconditioner\n");
643
644 if (opts->mat_vect_prod == DIRECT)
645 indPrecond_direct(sys, indsys, 2*PI*freq);
646 else if (opts->mat_vect_prod == MULTIPOLE)
647 indPrecond(sys, indsys, 2*PI*freq);
648 else {
649 fprintf(stderr, "Internal error: mat_vect_prod == %d\n",
650 opts->mat_vect_prod);
651 exit(1);
652 }
653 }
654 else if (indsys->precond_type == SPARSE) {
655 printf("Forming sparse matrix preconditioner..\n");
656 fill_spPre(sys, indsys, 2*PI*freq);
657
658 if (m == 0) {
659 if (indsys->opts->debug == ON)
660 printf("Number of nonzeros before factoring: %d\n",
661 spElementCount(indsys->sparMatrix));
662
663 /* Reorder and Factor the matrix */
664 err = spOrderAndFactor(indsys->sparMatrix, NULL, 1e-3, 0.0, 1);
665
666 if (indsys->opts->debug == ON)
667 printf("Number of fill-ins after factoring: %d\n",
668 spFillinCount(indsys->sparMatrix));
669 }
670 else
671 err = spFactor(indsys->sparMatrix);
672
673 if (err != 0) {
674 fprintf(stderr,"Error on factor: %d\n",err);
675 exit(1);
676 }
677 }
678 }
679 else if (opts->soln_technique == LUDECOMP) {
680 if (!dont_form_Z) {
681 printf("Performing LU Decomposition...\n");
682 cx_ludecomp(MtZM, num_mesh, FALSE);
683 printf("Done.\n");
684 }
685 else {
686 /* let's form a sparse version since fmin=0 and the L matrix = 0 */
687 create_sparMatrix(indsys);
688 printf("Filling sparse version M R Mt...");
689 fill_diagR(indsys);
690
691 /* dump matrix to disk if requested */
692 if (indsys->opts->dumpMats & MZMt) {
693 printf("saving sparse MZMt to MZMt.mat...\n");
694 concat4(outfname,"MZMt",indsys->opts->suffix,".mat");
695 if (spFileMatrix(indsys->sparMatrix, outfname, "MZMt", 0, 1, 1) == 0)
696 fprintf(stderr,"saving sparse matrix failed\n");
697 }
698
699 if (indsys->opts->debug == ON)
700 printf("Number of nonzeros before factoring: %d\n",
701 spElementCount(indsys->sparMatrix));
702
703 /* Reorder and Factor the matrix */
704 /* since w = 0, this matrix is real but all complex computations
705 are done since that is how the sparse package was compiled.
706 But changing it doesn't help that much. Must be dominated by
707 reordering and such.*/
708 printf("Factoring the sparse matrix...\n");
709 err = spOrderAndFactor(indsys->sparMatrix, NULL, 1e-3, 0.0, 1);
710
711 printf("Done factoring.\n");
712 if (indsys->opts->debug == ON)
713 printf("Number of fill-ins after factoring: %d\n",
714 spFillinCount(indsys->sparMatrix));
715 }
716 }
717 else {
718 fprintf(stderr, "Internal error: Unknown type of soln_technique: %d\n",
719 opts->soln_technique);
720 exit(1);
721 }
722
723 stoptimer;
724 ftimes[3] += dtime;
725
726 if (indsys->opts->debug == ON)
727 printf("Time spent on forming Precond: %lg\n",dtime);
728
729 starttimer;
730 for(ext = get_next_ext(indsys->externals), i=0; ext != NULL;
731 ext = get_next_ext(ext->next),i++) {
732 printf("conductor %d from node %s\n",i, get_a_name(ext->source));
733
734 /* initialize b */
735 for(j = 0; j < num_mesh; j++)
736 b[j] = x0[j] = CXZERO;
737
738 fill_b(ext, b);
739
740 sprintf(fname, "b%d_%d",m,i);
741 if (opts->dumpMats != OFF)
742 savecmplx(fb, fname, &b, 1, num_mesh);
743
744 maxiters = MIN(user_maxiters, num_mesh+2);
745
746 #if OPCNT == ON
747 maxiters = 2;
748 #endif
749
750 if (opts->soln_technique == ITERATIVE) {
751 printf("Calling gmres...\n");
752 if (opts->mat_vect_prod == MULTIPOLE)
753 gmres(MtZM, b, x0, inner, SetupComputePsi, num_mesh, maxiters, tol,
754 sys, chglist, 2*PI*freq, R, indsys, i);
755 else
756 gmres(MtZM, b, x0, inner, directmatvec, num_mesh, maxiters, tol,
757 sys, chglist, 2*PI*freq, R, indsys, i);
758
759 if (indsys->precond_type == LOC) {
760 multPrecond(indsys->Precond, x0, vect, num_mesh);
761 for (k=0; k < num_mesh; k++)
762 x0[k] = vect[k];
763 }
764 else if (indsys->precond_type == SPARSE)
765 spSolve(indsys->sparMatrix, x0, x0);
766
767 }
768 else
769 if (!dont_form_Z)
770 cx_lu_solve(MtZM, x0, b, num_mesh);
771 else
772 spSolve(indsys->sparMatrix, b, x0);
773
774 if (opts->dumpMats & GRIDS)
775 /* Do stuff to look at groundplane current distribution */
776 makegrids(indsys, x0, ext->Yindex, m);
777
778 /* Extract appropriate elements of x0 that correspond to final Yc */
779 extractYcol(indsys->FinalY, x0, ext, indsys->externals);
780
781 if (indsys->opts->debug == ON) {
782 /* SRW -- this is binary data */
783 fptemp = fopen("Ytemp.mat", "wb");
784 if (fptemp == NULL) {
785 printf("couldn't open file %s\n","Ytemp.mat");
786 exit(1);
787 }
788 strcpy(fname,"Ycond");
789 sprintf(tempstr, "%d", m);
790 strcat(fname,tempstr);
791 savecmplx(fptemp, fname, indsys->FinalY, num_extern, num_sub_extern);
792 fclose(fptemp);
793 }
794
795 }
796 stoptimer;
797 ftimes[4] += dtime;
798
799 if (i != indsys->num_sub_extern) {
800 fprintf(stderr, "Huh? columns calculated = %d and num extern = %d\n",
801 i, indsys->num_sub_extern);
802 exit(1);
803 }
804
805 if (indsys->opts->debug == ON) {
806 /* dump matrix to matlab file */
807 dump_to_Ycond(fp, m, indsys);
808 }
809
810 /* dump matrix to text file */
811 if (num_extern == num_sub_extern) {
812 fprintf(fp3, "Impedance matrix for frequency = %lg %d x %d\n ", freq,
813 num_extern, num_extern);
814 cx_invert(indsys->FinalY, num_extern);
815 }
816 else
817 fprintf(fp3, "ADMITTANCE matrix for frequency = %lg %d x %d\n ", freq,
818 num_extern, num_sub_extern);
819
820 cx_dumpMat_totextfile(fp3, indsys->FinalY, num_extern, num_sub_extern);
821 fflush(fp3);
822 }
823
824 if (indsys->opts->debug == ON)
825 fclose(fp);
826
827 fclose(fp3);
828 if (opts->dumpMats != OFF)
829 fclose(fb);
830
831 printf("\nAll impedance matrices dumped to file Zc%s.mat\n\n",opts->suffix);
832 totaltime = 0;
833 for(i = 0; i < TIMESIZE; i++)
834 totaltime += ftimes[i];
835
836 if (indsys->opts->debug == ON) {
837 printf("Calls to exact_mutual: %15d\n",num_exact_mutual);
838 printf(" fourfils: %15d\n",num_fourfil);
839 printf(" mutualfil: %15d\n",num_mutualfil);
840 printf("Number found in table: %15d\n",num_found);
841 printf("Number perpendicular: %15d\n",num_perp);
842 printf("\n");
843 }
844
845 printf("Times: Read geometry %lg\n",ftimes[0]);
846 printf(" Multipole setup %lg\n",ftimes[5]);
847 printf(" Scanning graph %lg\n",ftimes[6]);
848 printf(" Form A M and Z %lg\n",ftimes[1]);
849 printf(" form M'ZM %lg\n",ftimes[2]);
850 printf(" Form precond %lg\n",ftimes[3]);
851 printf(" GMRES time %lg\n",ftimes[4]);
852 printf(" Total: %lg\n",totaltime);
853
854 #ifdef MATTDEBUG
855 /* print memory bins */
856 for(i = 0; i < 1001; i++)
857 fprintf(stderr, "%d\n", membins[i]);
858 #endif
859
860 return (0);
861 }
862
863 /* this will divide a rectangular segment into many filaments */
assignFil(seg,num_fils,chgptr)864 charge *assignFil(seg, num_fils, chgptr)
865 SEGMENT *seg;
866 int *num_fils;
867 charge *chgptr;
868 {
869
870 int i,j,k;
871 double x,y,z, delw, delh;
872 double hx, hy, hz;
873 int temp, counter;
874 int Hinc, Winc;
875 double Hdiv, Wdiv;
876 double height, width;
877 double h_from_edge, w_from_edge, min_height, min_width;
878 FILAMENT *tempfil, *filptr;
879 int countfils;
880 double wx, wy, wz, mag; /* direction of width */
881 NODES *node0, *node1;
882 charge *clast, *chg;
883 charge *ctemp; /* temporary place so can allocate all the charges at once*/
884 surface *dummysurf;
885 int indices[4], row, col;
886 double r_width, r_height;
887 /* ratio of element sizes (for geometrically increasing size)*/
888
889 Hinc = seg->hinc;
890 Winc = seg->winc;
891 r_height = seg->r_height;
892 r_width = seg->r_width;
893
894 clast = chgptr;
895
896 seg->num_fils = Winc*Hinc;
897 seg->filaments = (FILAMENT *)MattAlloc(Hinc*Winc,sizeof(FILAMENT));
898
899 ctemp = (charge *)MattAlloc(Hinc*Winc, sizeof(charge));
900
901 /* clear pchg as a marker for checking later */
902 for(i = 0; i < Hinc*Winc; i++)
903 seg->filaments[i].pchg = NULL;
904
905 dummysurf = (surface *)MattAlloc(1, sizeof(surface));
906 dummysurf->type = CONDTR;
907 dummysurf->next = NULL;
908 dummysurf->prev = NULL;
909
910 /* To make the filaments have geometrically decreasing areas */
911 /* Hdiv and Wdiv are 1/(smallest Width) */
912
913 if (fabs(1.0 - r_height) < EPS)
914 Hdiv = Hinc;
915 else {
916 temp = Hinc/2;
917 Hdiv = 2*(1.0 - pow(r_height, (double)temp))/(1.0 - r_height);
918 if (Hinc%2 == 1) Hdiv += pow(r_height, (double)(temp));
919 }
920
921 if (fabs(1.0 - r_width) < EPS)
922 Wdiv = Winc;
923 else {
924 temp = Winc/2;
925 Wdiv = 2*(1.0 - pow(r_width, (double)temp))/(1.0 - r_width);
926 if (Winc%2 == 1) Wdiv += pow(r_width, (double)(temp));
927 }
928
929 node0 = seg->node[0];
930 node1 = seg->node[1];
931 /* determine direction of width */
932 if (seg->widthdir != NULL) {
933 wx = seg->widthdir[XX];
934 wy = seg->widthdir[YY];
935 wz = seg->widthdir[ZZ];
936 }
937 else {
938 /* default for width direction is in x-y plane perpendic to length*/
939 /* so do cross product with unit z*/
940 wx = -(node1->y - node0->y)*1.0;
941 wy = (node1->x - node0->x)*1.0;
942 wz = 0;
943 if ( fabs(wx/seg->length) < EPS && fabs(wy/seg->length) < EPS) {
944 /* if all of x-y is perpendic to length, then choose x direction */
945 wx = 1.0;
946 wy = 0;
947 }
948 mag = sqrt(wx*wx + wy*wy + wz*wz);
949 wx = wx/mag;
950 wy = wy/mag;
951 wz = wz/mag;
952 }
953 /* height direction perpendicular to both */
954 hx = -wy*(node1->z - node0->z) + (node1->y - node0->y)*wz;
955 hy = -wz*(node1->x - node0->x) + (node1->z - node0->z)*wx;
956 hz = -wx*(node1->y - node0->y) + (node1->x - node0->x)*wy;
957 mag = sqrt(hx*hx + hy*hy + hz*hz);
958 hx = hx/mag;
959 hy = hy/mag;
960 hz = hz/mag;
961
962 filptr = seg->filaments;
963 counter = 0;
964
965 /* this will fill the 'filament' array. It uses symmetry wrt z and y */
966 /* it generates the four corner fils first, then the next one in... */
967 /* 6/25/93 - added stuff to place fils in filament array so that adjacent */
968 /* fils are adjacent in the array. This will make the meshes */
969 /* in M consist of fils that are near each other. */
970 h_from_edge = 0.0;
971 min_height = 1.0/Hdiv; /* fil of smallest height */
972 for(i = 0; i < Hinc/2 || (Hinc%2 == 1 && i == Hinc/2); i++) {
973
974 /* height of the filaments for this row */
975 height = (seg->height/Hdiv)*pow(r_height, (double)i);
976 /* delh = (seg->height)*(0.5 - 1.0/Hdiv*
977 ( (1 - pow(r_height, (double)(i+1)))/(1-r_height)
978 - 0.5*pow(r_height, (double) i) )
979 );
980 */
981 if (i == 0)
982 h_from_edge += min_height/2;
983 else
984 h_from_edge += min_height/2*pow(r_height,(double)(i-1))*(1+r_height);
985
986 delh = (seg->height)*(0.5 - h_from_edge);
987
988 if (delh < 0.0 && fabs(delh/seg->height) > EPS)
989 printf("uh oh, delh < 0. delh/height = %lg\n", delh/seg->height);
990
991 w_from_edge = 0;
992 min_width = 1.0/Wdiv;
993 for(j = 0; j < Winc/2 || (Winc%2 == 1 && j == Winc/2); j++) {
994 width = (seg->width/Wdiv)*pow(r_width, (double)j );
995 /*delw = (seg->width)*
996 (0.5 - 1.0/Wdiv*
997 ( (1 - pow(r_width, (double)(j+1)))/(1 - r_width)
998 - 0.5*pow(r_width, (double) j) )
999 );
1000 */
1001
1002 if (j == 0)
1003 w_from_edge += min_width/2;
1004 else
1005 w_from_edge += min_width/2*pow(r_width,(double)(j-1))*(1+r_width);
1006
1007 delw = (seg->width)*(0.5 - w_from_edge);
1008
1009 if (delw < 0.0 && fabs(delw/seg->width) > EPS)
1010 printf("uh oh, delw < 0. delw/width = %lg\n", delw/seg->width);
1011 /* tempfil = filptr; */
1012 countfils = 0;
1013
1014 row = i;
1015 col = j;
1016 if (row%2 == 1)
1017 col = (Winc - 1) - col;
1018 indices[countfils] = col + Winc*row;
1019 filptr = &(seg->filaments[indices[countfils]]);
1020
1021 filptr->x[0] = node0->x + hx*delh + wx*delw;
1022 filptr->x[1] = node1->x + hx*delh + wx*delw;
1023 filptr->y[0] = node0->y + hy*delh + wy*delw;
1024 filptr->y[1] = node1->y + hy*delh + wy*delw;
1025 filptr->z[0] = node0->z + hz*delh + wz*delw;
1026 filptr->z[1] = node1->z + hz*delh + wz*delw;
1027 filptr->pchg = &ctemp[counter++];
1028 /* filptr++; */
1029 countfils++;
1030
1031 /* do symmetric element wrt y */
1032 if(j != Winc/2) {
1033 row = i;
1034 col = (Winc - 1) - j;
1035 if (row%2 == 1)
1036 col = (Winc - 1) - col;
1037 indices[countfils] = col + Winc*row;
1038 filptr = &(seg->filaments[indices[countfils]]);
1039
1040 filptr->x[0] = node0->x + hx*delh - wx*delw;
1041 filptr->x[1] = node1->x + hx*delh - wx*delw;
1042 filptr->y[0] = node0->y + hy*delh - wy*delw;
1043 filptr->y[1] = node1->y + hy*delh - wy*delw;
1044 filptr->z[0] = node0->z + hz*delh - wz*delw;
1045 filptr->z[1] = node1->z + hz*delh - wz*delw;
1046 filptr->pchg = &ctemp[counter++];
1047 /* filptr++; */
1048 countfils++;
1049 }
1050
1051 /* wrt z */
1052 if(i != Hinc/2) {
1053 row = (Hinc - 1) - i;
1054 col = j;
1055 if (row%2 == 1)
1056 col = (Winc - 1) - col;
1057 indices[countfils] = col + Winc*row;
1058 filptr = &(seg->filaments[indices[countfils]]);
1059
1060 filptr->x[0] = node0->x - hx*delh + wx*delw;
1061 filptr->x[1] = node1->x - hx*delh + wx*delw;
1062 filptr->y[0] = node0->y - hy*delh + wy*delw;
1063 filptr->y[1] = node1->y - hy*delh + wy*delw;
1064 filptr->z[0] = node0->z - hz*delh + wz*delw;
1065 filptr->z[1] = node1->z - hz*delh + wz*delw;
1066 filptr->pchg = &ctemp[counter++];
1067 filptr++;
1068 countfils++;
1069 }
1070
1071 /* wrt z and y */
1072 if( i != Hinc/2 && j != Winc/2) {
1073 row = (Hinc - 1) - i;
1074 col = (Winc - 1) - j;
1075 if (row%2 == 1)
1076 col = (Winc - 1) - col;
1077 indices[countfils] = col + Winc*row;
1078 filptr = &(seg->filaments[indices[countfils]]);
1079
1080 filptr->x[0] = node0->x - hx*delh - wx*delw;
1081 filptr->x[1] = node1->x - hx*delh - wx*delw;
1082 filptr->y[0] = node0->y - hy*delh - wy*delw;
1083 filptr->y[1] = node1->y - hy*delh - wy*delw;
1084 filptr->z[0] = node0->z - hz*delh - wz*delw;
1085 filptr->z[1] = node1->z - hz*delh - wz*delw;
1086 filptr->pchg = &ctemp[counter++];
1087 /* filptr++; */
1088 countfils++;
1089 }
1090
1091 for(k = 0; k < countfils; k++) {
1092 tempfil = &(seg->filaments[indices[k]]);
1093 tempfil->length = seg->length;
1094 tempfil->area = width*height;
1095 tempfil->width = width;
1096 tempfil->height = height;
1097 tempfil->filnumber = (*num_fils)++;
1098 tempfil->segm = seg;
1099
1100 tempfil->lenvect[XX] = tempfil->x[1] - tempfil->x[0];
1101 tempfil->lenvect[YY] = tempfil->y[1] - tempfil->y[0];
1102 tempfil->lenvect[ZZ] = tempfil->z[1] - tempfil->z[0];
1103 /* do stuff for multipole */
1104 /* make linked list entry */
1105 chg = tempfil->pchg;
1106 clast->next = chg;
1107 clast = chg;
1108 /* fill charge structure */
1109 chg->max_diag = chg->min_diag = tempfil->length;
1110 chg->x = (tempfil->x[0] + tempfil->x[1])/2.0;
1111 chg->y = (tempfil->y[0] + tempfil->y[1])/2.0;
1112 chg->z = (tempfil->z[0] + tempfil->z[1])/2.0;
1113 chg->surf = dummysurf;
1114 chg->dummy = FALSE;
1115 chg->fil = tempfil;
1116
1117 /* tempfil++; */
1118 }
1119
1120 }
1121 }
1122
1123 i = 0;
1124 while(i < Hinc*Winc && seg->filaments[i].pchg != NULL)
1125 i++;
1126
1127 if (i != Hinc*Winc) {
1128 fprintf(stderr, "Hey, not all filaments created in assignfil()! \n");
1129 exit(1);
1130 }
1131
1132 return clast;
1133 }
1134
1135
MatrixAlloc(rows,cols,size)1136 double **MatrixAlloc(rows, cols, size)
1137 int rows, cols, size;
1138 {
1139
1140 double **temp;
1141 int i;
1142
1143 temp = (double **)MattAlloc(rows,sizeof(double *));
1144 if (temp == NULL) {
1145 printf("not enough space for matrix allocation\n");
1146 exit(1);
1147 }
1148
1149 for(i = 0; i < rows; i++)
1150 temp[i] = (double *)MattAlloc(cols,size);
1151
1152 if (temp[rows - 1] == NULL) {
1153 printf("not enough space for matrix allocation\n");
1154 exit(1);
1155 }
1156 return(temp);
1157 }
1158
fillA(indsys)1159 void fillA(indsys)
1160 SYS *indsys;
1161 {
1162 SEGMENT *seg;
1163 NODES *node1, *node2, *node;
1164 MELEMENT **Alist;
1165 int i, counter;
1166 FILAMENT *fil;
1167
1168 indsys->Alist = (MELEMENT **)MattAlloc(indsys->num_real_nodes,
1169 sizeof(MELEMENT *));
1170
1171 pick_ground_nodes(indsys);
1172
1173 Alist = indsys->Alist;
1174 counter = 1; /* ground is chosen already */
1175
1176 for(seg = indsys->segment; seg != NULL; seg = seg->next) {
1177 node1 = getrealnode(seg->node[0]);
1178 node2 = getrealnode(seg->node[1]);
1179 if (node1->index == -1) {
1180 node1->index = counter++;
1181 Alist[node1->index] = NULL;
1182 }
1183 if (node2->index == -1) {
1184 node2->index = counter++;
1185 Alist[node2->index] = NULL;
1186 }
1187 for(i = 0; i < seg->num_fils; i++) {
1188 fil = &seg->filaments[i];
1189 Alist[node1->index] = insert_in_list(make_melement(fil->filnumber,
1190 fil, 1),
1191 Alist[node1->index]);
1192 Alist[node2->index] = insert_in_list(make_melement(fil->filnumber,
1193 fil, -1),
1194 Alist[node2->index]);
1195 }
1196 }
1197
1198 if (counter != indsys->num_real_nodes - indsys->num_trees + 1) {
1199 fprintf(stderr,"Internal error when forming A: counter %d != num_real_nodes %d\n",
1200 counter, indsys->num_real_nodes);
1201 }
1202
1203 }
1204
1205 /* this fills the kircoff's voltage law matrix (Mesh matrix) */
1206 /* it maps a matrix of mesh currents to branch currents */
1207 /* it might actually be what some think of as the transpose of M */
1208 /* Here, M*Im = Ib where Im are the mesh currents, and Ib the branch */
1209 /* 6/92 I added Mlist which is a vector of linked lists to meshes.
1210 This replaces M. But I keep M around for checking things in matlab. */
old_fillM(indsys)1211 void old_fillM(indsys)
1212 SYS *indsys;
1213 {
1214 }
1215
fillZ(indsys)1216 void fillZ(indsys)
1217 SYS *indsys;
1218 {
1219 int i, j, k, m;
1220 FILAMENT *fil_j, *fil_m;
1221 int filnum_j, filnum_m;
1222 double w;
1223 SEGMENT *seg1, *seg2;
1224 double **Z, *R, freq;
1225 int num_segs;
1226
1227 Z = indsys->Z;
1228 R = indsys->R;
1229
1230 for(seg1 = indsys->segment; seg1 != NULL; seg1 = seg1->next) {
1231 for(j = 0; j < seg1->num_fils; j++) {
1232 fil_j = &(seg1->filaments[j]);
1233 filnum_j = fil_j->filnumber;
1234 #if SUPERCON == ON
1235 R[filnum_j] = fil_j->length*seg1->r1/fil_j->area;
1236 #else
1237 R[filnum_j] = resistance(fil_j, seg1->sigma);
1238 #endif
1239 if (indsys->opts->mat_vect_prod != MULTIPOLE
1240 && !indsys->dont_form_Z) {
1241 for(seg2 = indsys->segment; seg2 != NULL; seg2 = seg2->next) {
1242 for(m = 0; m < seg2->num_fils; m++) {
1243 fil_m = &(seg2->filaments[m]);
1244 filnum_m = fil_m->filnumber;
1245 if (filnum_m == filnum_j) {
1246 Z[filnum_m][filnum_m] = selfterm(fil_m); /* do self-inductance */
1247 #if SUPERCON == ON
1248 if (seg1->lambda != 0.0)
1249 Z[filnum_m][filnum_m] += seg1->r2*fil_j->length/fil_j->area;
1250 #endif
1251 }
1252 else
1253 if (filnum_m > filnum_j) /*we haven't done it yet */
1254
1255 Z[filnum_m][filnum_j]
1256 = Z[filnum_j][filnum_m] = mutual(fil_j, fil_m);
1257 }
1258 }
1259 } /* end if (multipole or dont_form_Z) */
1260 }
1261 }
1262 }
1263
1264 #if SUPERCON == ON
1265 /* Have to reset the diag elements for each omega, if superconductor
1266 * (lambda != 0) and sigma != 0.
1267 */
fillZ_diag(SYS * indsys,double omega)1268 void fillZ_diag(SYS *indsys, double omega)
1269 {
1270 int i, j;
1271 FILAMENT *fil_j;
1272 int filnum_j;
1273 SEGMENT *seg1;
1274 double **Z, *R;
1275 double tmp, tmp1, dnom, r1, r2;
1276
1277 Z = indsys->Z;
1278 R = indsys->R;
1279
1280 for(seg1 = indsys->segment; seg1 != NULL; seg1 = seg1->next) {
1281 if (seg1->lambda != 0.0 && seg1->sigma != 0.0) {
1282 /* segment is a superconductor with frequency dependent terms */
1283 tmp = MU0*seg1->lambda*seg1->lambda;
1284 tmp1 = tmp*omega;
1285 dnom = tmp1*seg1->sigma;
1286 dnom = dnom*dnom + 1.0;
1287 seg1->r1 = r1 = seg1->sigma*tmp1*tmp1/dnom;
1288 seg1->r2 = r2 = tmp/dnom;
1289 for(j = 0; j < seg1->num_fils; j++) {
1290 fil_j = &(seg1->filaments[j]);
1291 filnum_j = fil_j->filnumber;
1292 R[filnum_j] = fil_j->length*r1/fil_j->area;
1293 Z[filnum_j][filnum_j] = selfterm(fil_j) +
1294 r2*fil_j->length/fil_j->area;
1295 }
1296 }
1297 }
1298 }
1299
1300 /*
1301 * Theory:
1302 * In normal metal: (1) del X H = -i*omega*mu*sigma * H
1303 * In superconductor: (2) del X H = (1/lambda)^2 * H
1304 *
1305 * In fasthenry, (1) is solved, so the game is to replace sigma in (1)
1306 * with a complex variable that includes and reduces to (2). We choose
1307 *
1308 * sigma_prime = sigma + i/(omega*mu*lambda^2)
1309 *
1310 * Then, using sigma_prime in (1) rather than sigma, one obtains an
1311 * expression that reduces to (2) at omega = 0, yet retains properties
1312 * of (1). This is the two-fluid model, where the sigma in sigma_prime
1313 * represents the conductivity due to unpaired electrons.
1314 *
1315 * Since sigma_prime blows up at omega = 0, we work with the impedance,
1316 * which we take as z = r1 + i*omega*r2 = i/sigma_prime. The r1 and
1317 * r2 variables are thus
1318 *
1319 * (3) r1 = sigma*(omega*mu*lambda^2)^2
1320 * --------------------------------
1321 * (sigma*omega*mu*lambda^2)^2 + 1
1322 *
1323 * (4) r2 = mu*lambda^2
1324 * --------------------------------
1325 * (sigma*omega*mu*lambda^2)^2 + 1
1326 */
1327
1328 /* Set the r1, r2 entries to the appropriate values. The r1 entry
1329 * comes from the real value of sigma (1/sigma for normal metals).
1330 * The r2 entry comes from the imaginary part of sigma, which reduces
1331 * to MU0*lambda^2 when the real part of sigma (default 0 for
1332 * superconductors) is negligible, and is 0 for normal metals.
1333 * The r2 term contributes to the inductance matrix diagonal.
1334 */
set_rvals(SYS * indsys,double omega)1335 void set_rvals(SYS *indsys, double omega)
1336 {
1337 SEGMENT *seg1;
1338 double tmp, tmp1, dnom, r1, r2;
1339
1340 for(seg1 = indsys->segment; seg1 != NULL; seg1 = seg1->next) {
1341 if (seg1->lambda != 0.0) {
1342 /* segment is a superconductor */
1343 tmp = MU0*seg1->lambda*seg1->lambda;
1344 if (seg1->sigma != 0.0) {
1345 /* the terms become frequency dependent */
1346 tmp1 = tmp*omega;
1347 dnom = tmp1*seg1->sigma;
1348 dnom = dnom*dnom + 1.0;
1349 r1 = seg1->sigma*tmp1*tmp1/dnom;
1350 r2 = tmp/dnom;
1351 }
1352 else {
1353 r1 = 0.0;
1354 r2 = tmp;
1355 }
1356 }
1357 else {
1358 r1 = 1/seg1->sigma;
1359 r2 = 0.0;
1360 }
1361 seg1->r1 = r1;
1362 seg1->r2 = r2;
1363 }
1364 }
1365 #endif
1366
1367 /* calculates resistance of filament */
resistance(fil,sigma)1368 double resistance(fil, sigma)
1369 FILAMENT *fil;
1370 double sigma; /* conductivitiy */
1371 {
1372 return fil->length/(sigma * fil->area);
1373 }
1374
1375 /* mutual inductance functions moved to mutual.c */
1376
matherr(exc)1377 int matherr(exc)
1378 struct exception *exc;
1379 {
1380 printf("Err in math\n");
1381 return(0);
1382 }
1383
1384 /* This counts the nonblank lines of the file fp (unused) */
countlines(fp)1385 int countlines(fp)
1386 FILE *fp;
1387 {
1388 int count;
1389 char temp[MAXCHARS], *returnstring;
1390
1391 count = 0;
1392 while( fgets(temp,MAXCHARS, fp) != NULL)
1393 if ( notblankline(temp) ) count++;
1394
1395 return count;
1396 }
1397
1398 /* returns 1 if string contains a nonspace character */
notblankline(string)1399 static int notblankline(string)
1400 char *string;
1401 {
1402 while( *string!='\0' && isspace(*string))
1403 string++;
1404
1405 if (*string == '\0') return 0;
1406 else return 1;
1407 }
1408
1409 /* This saves various matrices to files and optionally calls fillA() if
1410 the incidence matrix, A, is requested */
savemats(indsys)1411 savemats(indsys)
1412 SYS *indsys;
1413 {
1414 int i, j;
1415 FILE *fp, *fp2;
1416 FILAMENT *tmpf;
1417 SEGMENT *tmps;
1418 double *dumb, *dumbx, *dumby, *dumbz;
1419 int num_mesh, num_fils, num_real_nodes;
1420 double *Mrow; /* one row of M */
1421 double **Z, *R;
1422 int machine, counter, nnz, nnz0;
1423 ind_opts *opts;
1424 MELEMENT *mesh;
1425 MELEMENT **Mlist = indsys->Mlist;
1426 MELEMENT **Alist;
1427
1428 num_mesh = indsys->num_mesh;
1429 num_fils = indsys->num_fils;
1430 num_real_nodes = indsys->num_real_nodes;
1431
1432 Mrow = (double *)MattAlloc(num_fils,sizeof(double));
1433 R = indsys->R;
1434 Z = indsys->Z;
1435 opts = indsys->opts;
1436
1437 if (opts->dumpMats & DUMP_M) {
1438 /* count non-zero entries */
1439 nnz = nnz_inM(indsys->Mlist, num_mesh);
1440
1441 if (opts->kind & TEXT) {
1442 concat4(outfname,"M",opts->suffix,".dat");
1443 /* SRW -- this is ascii data */
1444 if ( (fp2 = fopen(outfname,"w")) == NULL) {
1445 printf("Couldn't open file\n");
1446 exit(1);
1447 }
1448 dump_M_to_text(fp2, Mlist, num_mesh, nnz);
1449 fclose(fp2);
1450 }
1451
1452 if (opts->kind & MATLAB) {
1453 concat4(outfname,"M",opts->suffix,".mat");
1454 /* SRW -- this is binary data */
1455 if ( (fp = fopen(outfname,"wb")) == NULL) {
1456 printf("Couldn't open file\n");
1457 exit(1);
1458 }
1459 dump_M_to_matlab(fp, Mlist, num_mesh, nnz, "M");
1460 fclose(fp);
1461 }
1462 }
1463
1464 if (opts->dumpMats & DUMP_A) {
1465
1466 /* fill the A matrix */
1467 fillA(indsys);
1468 Alist = indsys->Alist;
1469
1470 /* count non-zero entries */
1471 nnz = nnz_inM(indsys->Alist, num_real_nodes);
1472
1473 /* how many non-zeros in ground node 0? */
1474 nnz0 = nnz_inM(indsys->Alist, 1);
1475
1476 /* now dump it to a file without ground node */
1477
1478 if (opts->kind & TEXT) {
1479 concat4(outfname,"A",opts->suffix,".dat");
1480 /* SRW -- this is ascii data */
1481 if ( (fp2 = fopen(outfname,"w")) == NULL) {
1482 printf("Couldn't open file\n");
1483 exit(1);
1484 }
1485 dump_M_to_text(fp2, &Alist[1], num_real_nodes - 1, nnz - nnz0);
1486 fclose(fp2);
1487 }
1488
1489 if (opts->kind & MATLAB) {
1490 concat4(outfname,"A",opts->suffix,".mat");
1491 /* SRW -- this is binary data */
1492 if ( (fp = fopen(outfname,"wb")) == NULL) {
1493 printf("Couldn't open file\n");
1494 exit(1);
1495 }
1496 dump_M_to_matlab(fp, &Alist[1], num_real_nodes - 1, nnz - nnz0,"A");
1497 fclose(fp);
1498 }
1499 }
1500
1501 /* save text files */
1502 if (opts->kind & TEXT && opts->dumpMats & DUMP_RL) {
1503 /* SRW -- this is ascii data */
1504 /*
1505 if ( (fp2 = fopen("M.dat","w")) == NULL) {
1506 printf("Couldn't open file\n");
1507 exit(1);
1508 }
1509
1510 for(i = 0; i < num_fils; i++)
1511 Mrow[i] = 0;
1512
1513 for(i = 0; i < num_mesh; i++) {
1514 fillMrow(indsys->Mlist, i, Mrow);
1515 dumpVec_totextfile(fp2, Mrow, num_fils);
1516 }
1517
1518 fclose(fp2);
1519 */
1520
1521 concat4(outfname,"R",opts->suffix,".dat");
1522 /* SRW -- this is ascii data */
1523 if ( (fp2 = fopen(outfname,"w")) == NULL) {
1524 printf("Couldn't open file\n");
1525 exit(1);
1526 }
1527
1528 dumpVec_totextfile(fp2, R, num_fils);
1529
1530 fclose(fp2);
1531
1532 if (!indsys->dont_form_Z && indsys->opts->mat_vect_prod == DIRECT) {
1533 /* do imaginary part of Z */
1534
1535 concat4(outfname,"L",opts->suffix,".dat");
1536 /* SRW -- this is ascii data */
1537 if ( (fp2 = fopen(outfname,"w")) == NULL) {
1538 printf("Couldn't open file\n");
1539 exit(1);
1540 }
1541
1542 dumpMat_totextfile(fp2, Z, num_fils, num_fils);
1543
1544 fclose(fp2);
1545
1546 }
1547 else {
1548 if (indsys->dont_form_Z)
1549 printf("L matrix not dumped since L = 0 since fmin = 0\n");
1550 else
1551 printf("L matrix not dumped. Run with mat_vect_prod = DIRECT if this is desired\n");
1552 }
1553 }
1554
1555 /* Save Matlab files */
1556
1557 if (opts->kind & MATLAB && opts->dumpMats & DUMP_RL) {
1558 concat4(outfname,"RL",opts->suffix,".mat");
1559 /* SRW -- this is binary data */
1560 if ( (fp = fopen(outfname,"wb")) == NULL) {
1561 printf("Couldn't open file\n");
1562 exit(1);
1563 }
1564
1565 machine = 1000;
1566 #ifdef DEC
1567 machine = 0000;
1568 #endif
1569
1570 dumbx = (double *)malloc(num_fils*sizeof(double));
1571 dumby = (double *)malloc(num_fils*sizeof(double));
1572 dumbz = (double *)malloc(num_fils*sizeof(double));
1573 dumb = (double *)malloc(num_fils*sizeof(double));
1574 if (dumb == NULL) {
1575 printf("no space for R\n");
1576 exit(1);
1577 }
1578
1579 /* save and print matrices */
1580
1581 /*
1582 for(i = 0; i < num_fils; i++)
1583 Mrow[i] = 0;
1584
1585 for(i = 0; i < num_mesh; i++) {
1586 fillMrow(indsys->Mlist, i, Mrow);
1587 savemat_mod(fp, machine+100, "M", num_mesh, num_fils, 0, Mrow,
1588 (double *)NULL, i, num_fils);
1589 }
1590 */
1591
1592 /* this only saves the real part (savemat_mod.c modified) */
1593 savemat_mod(fp, machine+100, "R", 1, num_fils, 0, R, (double *)NULL,
1594 0, num_fils);
1595
1596 if (!indsys->dont_form_Z && indsys->opts->mat_vect_prod == DIRECT) {
1597 /* do imaginary part of Z */
1598 for(i = 0; i < num_fils; i++) {
1599 savemat_mod(fp, machine+100, "L", num_fils, num_fils, 0, Z[i],
1600 (double *)NULL, i, num_fils);
1601 }
1602 }
1603 else {
1604 if (indsys->dont_form_Z)
1605 printf("L matrix not dumped since L = 0 since fmin = 0\n");
1606 else
1607 printf("L matrix not dumped. Run with mat_vect_prod = DIRECT if this is desired\n");
1608 }
1609
1610 /* save vector of filament areas */
1611 for(tmps = indsys->segment; tmps != NULL; tmps = tmps->next) {
1612 for(j = 0; j < tmps->num_fils; j++) {
1613 tmpf = &(tmps->filaments[j]);
1614 dumb[tmpf->filnumber] = tmpf->area;
1615 dumbx[tmpf->filnumber] = tmpf->x[0];
1616 dumby[tmpf->filnumber] = tmpf->y[0];
1617 dumbz[tmpf->filnumber] = tmpf->z[0];
1618 }
1619 }
1620 savemat_mod(fp, machine+0, "areas", num_fils, 1, 0, dumb, (double *)NULL,0,
1621 num_fils);
1622 savemat_mod(fp, machine+0, "pos", num_fils, 3, 0, dumbx, (double *)NULL, 0,
1623 num_fils);
1624 savemat_mod(fp, machine+0, "pos", num_fils, 3, 0, dumby, (double *)NULL, 1,
1625 num_fils);
1626 savemat_mod(fp, machine+0, "pos", num_fils, 3, 0, dumbz, (double *)NULL, 1,
1627 num_fils);
1628
1629 free(dumb);
1630 free(dumbx);
1631 free(dumby);
1632 free(dumbz);
1633
1634 fclose(fp);
1635 }
1636 }
1637
savecmplx(fp,name,Z,rows,cols)1638 savecmplx(fp, name, Z, rows, cols)
1639 FILE *fp;
1640 char *name;
1641 CX **Z;
1642 int rows, cols;
1643 {
1644 int i,j;
1645 int machine;
1646
1647 machine = 1000;
1648 #ifdef DEC
1649 machine = 0000;
1650 #endif
1651
1652 /* this only saves the real part (savemat_mod.c modified) one byte per call*/
1653 for(i = 0; i < rows; i++)
1654 for(j = 0; j < cols; j++)
1655 savemat_mod(fp, machine+100, name, rows, cols, 1, &Z[i][j].real,
1656 (double *)NULL, i+j, 1);
1657
1658 /* do imaginary part of Z */
1659 for(i = 0; i < rows; i++)
1660 for(j = 0; j < cols; j++)
1661 savemat_mod(fp, machine+100, name, rows, cols, 0, &Z[i][j].imag,
1662 (double *)NULL, 1, 1);
1663 }
1664
1665 /* saves a complex matrix more efficiently? */
savecmplx2(fp,name,Z,rows,cols)1666 savecmplx2(fp, name, Z, rows, cols)
1667 FILE *fp;
1668 char *name;
1669 CX **Z;
1670 int rows, cols;
1671 {
1672 int i,j;
1673 int machine;
1674 static double *temp;
1675 static int colmax = 0;
1676
1677 if (colmax < cols) {
1678 temp = (double *)malloc(cols*sizeof(double));
1679 colmax = cols;
1680 }
1681
1682 machine = 1000;
1683 #ifdef DEC
1684 machine = 0000;
1685 #endif
1686
1687 /* this only saves the real part (savemat_mod.c modified) one byte per call*/
1688 for(i = 0; i < rows; i++) {
1689 for(j = 0; j < cols; j++)
1690 temp[j] = Z[i][j].real;
1691 savemat_mod(fp, machine+100, name, rows, cols, 1, temp,
1692 (double *)NULL, i, cols);
1693 }
1694
1695 /* do imaginary part of Z */
1696 for(i = 0; i < rows; i++) {
1697 for(j = 0; j < cols; j++)
1698 temp[j] = Z[i][j].imag;
1699 savemat_mod(fp, machine+100, name, rows, cols, 0, temp,
1700 (double *)NULL, 1, cols);
1701 }
1702 }
1703
1704 /* This computes the product M*Z*Mt in a better way than oldformMZMt */
formMZMt(indsys)1705 formMZMt(indsys)
1706 SYS *indsys;
1707 {
1708 int m,n,p;
1709 double tempsum, tempR, tempsumR;
1710 static double *tcol = NULL; /* temporary storage for extra rows */
1711 int i,j, k, mesh, mesh2, nodeindx;
1712 int nfils, nmesh;
1713 MELEMENT *melem, *melem2;
1714 MELEMENT *mt, *mt2; /* elements of M transpose */
1715 double **M, **L, *R;
1716 CX **Zm, *tempZ;
1717 int rows, cols, num_mesh;
1718 MELEMENT **Mlist, **Mt;
1719
1720 Zm = indsys->MtZM;
1721 L = indsys->Z;
1722 R = indsys->R;
1723 M = indsys->M;
1724 nfils = rows = indsys->num_fils;
1725 nmesh = cols = num_mesh = indsys->num_mesh;
1726 Mlist = indsys->Mlist;
1727 Mt = indsys->Mtrans;
1728
1729 if (nmesh > nfils) {
1730 fprintf(stderr, "Uh oh, more meshes than filaments, I'm confused\n");
1731 exit(1);
1732 }
1733
1734 /* allocate a temporary column for Z*Mt */
1735 if (tcol == NULL)
1736 tcol = (double *)MattAlloc(rows, sizeof(double));
1737
1738 /* this does L*(Mt)_i where (Mt)_i is a single column (row) of Mt (M) and
1739 saves it in the temp space, tcol */
1740 for(mesh = 0; mesh < num_mesh; mesh++) {
1741 for(j = 0; j< nfils; j++)
1742 tcol[j] = 0;
1743 /* note, this next set of nested loops could be reversed, but I think
1744 this might be more efficient for pipelining, ? */
1745 for(melem = Mlist[mesh]; melem != NULL; melem = melem->mnext)
1746 for(j = 0; j < nfils; j++)
1747 tcol[j] += L[j][melem->filindex]*melem->sign;
1748 for(mesh2 = 0; mesh2 < num_mesh; mesh2++) {
1749 tempsum = 0;
1750 for(melem2 = Mlist[mesh2]; melem2 != NULL; melem2 = melem2->mnext)
1751 tempsum += melem2->sign*tcol[melem2->filindex];
1752 Zm[mesh2][mesh].imag = tempsum;
1753 }
1754 }
1755
1756 /* R is diagonal, so M*R*Mt is easy */
1757 for(i = 0; i < num_mesh; i++)
1758 for(j = 0; j < num_mesh; j++)
1759 Zm[i][j].real = 0;
1760
1761 for(i = 0; i < nfils; i++)
1762 for(mt = Mt[i]; mt != NULL; mt = mt->mnext) {
1763 tempR = R[i]*mt->sign;
1764 for(mt2 = Mt[i]; mt2 != NULL; mt2 = mt2->mnext)
1765 Zm[mt2->filindex][mt->filindex].real += mt2->sign*tempR;
1766 }
1767 }
1768
oldformMZMt(indsys)1769 oldformMZMt(indsys)
1770 SYS *indsys;
1771 {
1772 int m,n,p;
1773 double tempsum;
1774 static CX **trows = NULL; /* temporary storage for extra rows */
1775 int i,j, k, mesh, nodeindx;
1776 int nfils, nmesh;
1777 MELEMENT *melem, *melem2;
1778 double **M, **L, *R;
1779 CX **Zm;
1780 int rows, cols, num_mesh;
1781 MELEMENT **Mlist;
1782
1783 Zm = indsys->MtZM;
1784 L = indsys->Z;
1785 R = indsys->R;
1786 M = indsys->M;
1787 nfils = rows = indsys->num_fils;
1788 nmesh = cols = num_mesh = indsys->num_mesh;
1789 Mlist = indsys->Mlist;
1790
1791 if (nmesh > nfils) {
1792 fprintf(stderr, "Uh oh, more meshes than filaments, I'm confused\n");
1793 exit(1);
1794 }
1795
1796 /* allocate some extra rows so we can use Zm as temp space */
1797 if (rows > cols && trows == NULL)
1798 trows = (CX **)MatrixAlloc(rows - cols, cols, sizeof(CX));
1799
1800 /* this does L*Mt and saves it in Zm.real. This could be done better i
1801 think (use the fact that MZMt is symmetric)
1802 Also, could only track through each mesh list once, and store
1803 temporary values in Zm.real as we go along. (someday) */
1804 for(mesh = 0; mesh < num_mesh; mesh++) {
1805 for(j = 0; j < nfils; j++) {
1806 tempsum = 0;
1807 for(melem = Mlist[mesh]; melem != NULL; melem = melem->mnext)
1808 tempsum += L[j][melem->filindex]*melem->sign;
1809 if (j < nmesh)
1810 Zm[j][mesh].real = tempsum;
1811 else
1812 trows[j - nmesh][mesh].real = tempsum;
1813 }
1814 }
1815
1816 /* savecmplx(fp, "step1", Zm, num_mesh, num_mesh); */
1817
1818 /* this does M*(Zm.real) where Zm.real is hopefully L*Mt */
1819 for(mesh = 0; mesh < num_mesh; mesh++) {
1820 for(j = 0; j < nmesh; j++) {
1821 tempsum = 0;
1822 for(melem = Mlist[mesh]; melem != NULL; melem = melem->mnext) {
1823 if (melem->filindex < nmesh)
1824 tempsum += Zm[melem->filindex][j].real*melem->sign;
1825 else
1826 tempsum += trows[melem->filindex - nmesh][j].real*melem->sign;
1827 }
1828 Zm[mesh][j].imag = tempsum;
1829 }
1830 }
1831 /* savecmplx(fp, "step2", Zm, num_mesh, num_mesh); */
1832
1833 /* R is diagonal, so M*R*Mt is easy */
1834 for(i = 0; i < num_mesh; i++) {
1835 for(j = i; j < num_mesh; j++) {
1836 tempsum = 0;
1837 /* brute force search for elements */
1838 for(melem = Mlist[i]; melem != NULL; melem = melem->mnext) {
1839 for(melem2 = Mlist[j]; melem2 != NULL; melem2 = melem2->mnext) {
1840 if (melem2->filindex == melem->filindex)
1841 tempsum += melem->sign*R[melem->filindex]*melem2->sign;
1842 }
1843 }
1844 Zm[i][j].real = Zm[j][i].real = tempsum;
1845 }
1846 }
1847 /* don't free the temp space */
1848 /*
1849 if (rows > cols) {
1850 for(i = 0; i < rows - cols; i++)
1851 free(trows[i]);
1852 free(trows);
1853 }
1854 */
1855 }
1856
MattAlloc(number,size)1857 char* MattAlloc(number, size)
1858 int number, size;
1859 {
1860
1861 char *blah;
1862
1863 /* blah = (char *)malloc(number*size); */
1864 CALLOC(blah, number*size, char, OFF, IND);
1865
1866 if (blah == NULL) {
1867 fprintf(stderr, "MattAlloc: Couldn't get space. Needed %d\n",number*size);
1868 exit(1);
1869 }
1870
1871 return blah;
1872 }
1873
1874 /* forms the transpose of M. Makes a linked list of each row. Mlist is
1875 a linked list of its rows. */
1876 /* Note: This uses the same struct as Mlist but in reality, each linked list
1877 is composed of mesh indices, not fil indices. (filindex is a mesh index) */
formMtrans(indsys)1878 formMtrans(indsys)
1879 SYS *indsys;
1880 {
1881 int i, j, count;
1882 MELEMENT *m, *mt, *mt2, mtdum;
1883 MELEMENT **Mlist, **Mtrans, **Mtrans_check;
1884 int meshes, fils;
1885 int last_filindex;
1886 MELEMENT *create_melem();
1887
1888 fils = indsys->num_fils;
1889 meshes = indsys->num_mesh;
1890 Mlist = indsys->Mlist;
1891 Mtrans = indsys->Mtrans;
1892
1893 /* clear it (should be garbage only) */
1894 for(i = 0; i < fils; i++)
1895 Mtrans[i] = NULL;
1896
1897 for(j = 0; j < meshes; j++) {
1898 last_filindex = -1;
1899 for(m = Mlist[j]; m != NULL; m = m->mnext) {
1900 mt = make_melement(j, NULL, m->sign);
1901 Mtrans[m->filindex] = insert_in_list(mt,Mtrans[m->filindex]);
1902 if (last_filindex == m->filindex)
1903 printf("Internal Warning: Mesh %d contains the same filament multiple times\n",j);
1904 last_filindex = m->filindex;
1905 }
1906 }
1907
1908 #if 1==0
1909 /* this was the old inefficient way to from Mt */
1910
1911 /* allocated for comparison */
1912 Mtrans_check = (MELEMENT **)MattAlloc(fils,sizeof(MELEMENT *));
1913
1914 for(i = 0; i < fils; i++) {
1915 mt = &mtdum;
1916 /* scan through mesh list j looking for a filament i. (j,i) in M */
1917 for(j = 0; j < meshes; j++) {
1918 for(m = Mlist[j]; m->filindex < i && m->mnext != NULL; m = m->mnext)
1919 ;
1920 if (m->filindex == i) {
1921 mt->mnext = make_melement(j, NULL, m->sign);
1922 mt = mt->mnext;
1923 if (m->mnext != NULL && m->mnext->filindex == i) {
1924 for(count = 1; m->mnext->filindex == i; count++) {
1925 m = m->mnext;
1926 mt->mnext = make_melement(j, NULL, m->sign);
1927 mt = mt->mnext;
1928 }
1929 printf("Internal Warning: Mesh %d contains the same filament %d times\n",j,count);
1930 }
1931 }
1932 }
1933 mt->mnext = NULL;
1934 Mtrans_check[i] = mtdum.mnext;
1935 }
1936
1937 /* check old and new ways */
1938 for(i = 0; i < fils; i++) {
1939 for(mt = Mtrans[i], mt2 = Mtrans_check[i]; mt != NULL && mt2 != NULL;
1940 mt = mt->mnext, mt2 = mt2->mnext)
1941 if (mt->filindex != mt2->filindex || mt->sign != mt2->sign)
1942 printf("different: %d %d %d\n",i,mt->filindex, mt2->filindex);
1943 if (mt != mt2)
1944 printf("both not NULL: mt: %u mt2: %u\n", mt, mt2);
1945 }
1946 #endif
1947
1948 }
1949
compare_meshes(mesh1,mesh2)1950 compare_meshes(mesh1, mesh2)
1951 MELEMENT *mesh1, *mesh2;
1952 {
1953
1954 while(mesh1 != NULL && mesh2 != NULL && mesh1->filindex == mesh2->filindex && mesh1->sign == mesh2->sign) {
1955 mesh1 = mesh1->mnext;
1956 mesh2 = mesh2->mnext;
1957 }
1958
1959 if (mesh1 != NULL || mesh2 != NULL) {
1960 fprintf(stderr, "meshes don't match!\n");
1961 exit(1);
1962 }
1963 }
1964
cx_dumpMat_totextfile(fp,Z,rows,cols)1965 cx_dumpMat_totextfile(fp, Z, rows, cols)
1966 FILE *fp;
1967 CX **Z;
1968 int rows, cols;
1969 {
1970 int i, j;
1971
1972 for(i = 0; i < rows; i++) {
1973 for(j = 0; j < cols; j++)
1974 fprintf(fp, "%13.6lg %+13.6lgj ", Z[i][j].real, Z[i][j].imag);
1975 fprintf(fp, "\n");
1976 }
1977 return 0;
1978 }
1979
dumpMat_totextfile(fp,A,rows,cols)1980 dumpMat_totextfile(fp, A, rows, cols)
1981 FILE *fp;
1982 double **A;
1983 int rows, cols;
1984 {
1985 int i, j;
1986
1987 for(i = 0; i < rows; i++) {
1988 for(j = 0; j < cols; j++)
1989 fprintf(fp, "%13.6lg ", A[i][j]);
1990 fprintf(fp, "\n");
1991 }
1992 return 0;
1993 }
1994
dumpVec_totextfile(fp2,Vec,size)1995 dumpVec_totextfile(fp2, Vec, size)
1996 FILE *fp2;
1997 double *Vec;
1998 int size;
1999 {
2000 dumpMat_totextfile(fp2, &Vec, 1, size);
2001 }
2002
fillMrow(Mlist,mesh,Mrow)2003 fillMrow(Mlist, mesh, Mrow)
2004 MELEMENT **Mlist;
2005 int mesh;
2006 double *Mrow;
2007 {
2008 int i;
2009 MELEMENT *melem;
2010
2011 if (mesh != 0)
2012 /* remove non-zeros from last call */
2013 for(melem = Mlist[mesh - 1]; melem != NULL; melem=melem->mnext)
2014 Mrow[melem->filindex] = 0;
2015
2016 for(melem = Mlist[mesh]; melem != NULL; melem = melem->mnext)
2017 Mrow[melem->filindex] = melem->sign;
2018 }
2019
dump_to_Ycond(fp,cond,indsys)2020 dump_to_Ycond(fp, cond, indsys)
2021 FILE *fp;
2022 int cond;
2023 SYS *indsys;
2024 {
2025 static char fname[20], tempstr[5];
2026 int maxiters = indsys->opts->maxiters;
2027
2028 sprintf(tempstr, "%d", cond);
2029
2030 strcpy(fname,"Ycond");
2031 strcat(fname,tempstr);
2032 savecmplx(fp, fname, indsys->FinalY, indsys->num_extern,
2033 indsys->num_sub_extern);
2034
2035 strcpy(fname,"resids");
2036 strcat(fname,tempstr);
2037 saveCarray(fp, fname, indsys->resids, indsys->num_sub_extern, maxiters);
2038
2039 strcpy(fname,"resid_real");
2040 strcat(fname,tempstr);
2041 saveCarray(fp, fname, indsys->resid_real, indsys->num_sub_extern, maxiters);
2042
2043 strcpy(fname,"resid_imag");
2044 strcat(fname,tempstr);
2045 saveCarray(fp, fname, indsys->resid_imag, indsys->num_sub_extern, maxiters);
2046
2047 strcpy(fname,"niters");
2048 strcat(fname,tempstr);
2049 saveCarray(fp, fname, &(indsys->niters), 1, indsys->num_sub_extern);
2050
2051 }
2052
saveCarray(fp,fname,Arr,rows,cols)2053 saveCarray(fp, fname, Arr, rows, cols)
2054 FILE *fp;
2055 char *fname;
2056 double **Arr;
2057 int rows, cols;
2058 {
2059 int i;
2060 int machine;
2061
2062 machine = 1000;
2063 #ifdef DEC
2064 machine = 0000;
2065 #endif
2066
2067 for(i = 0; i < rows; i++) {
2068 savemat_mod(fp, machine+100, fname, rows, cols, 0, Arr[i], (double *)NULL,
2069 i, cols);
2070 }
2071 }
2072
nnz_inM(Mlist,num_mesh)2073 int nnz_inM(Mlist, num_mesh)
2074 MELEMENT **Mlist;
2075 int num_mesh;
2076 {
2077 int counter, i;
2078 MELEMENT *mesh;
2079
2080 counter = 0;
2081
2082 for(i = 0; i < num_mesh; i++)
2083 for(mesh = Mlist[i]; mesh != NULL; mesh = mesh->mnext)
2084 counter++;
2085
2086 return counter;
2087 }
2088
dump_M_to_text(fp,Mlist,num_mesh,nnz)2089 dump_M_to_text(fp, Mlist, num_mesh, nnz)
2090 FILE *fp;
2091 MELEMENT **Mlist;
2092 int num_mesh, nnz;
2093 {
2094 int counter, i;
2095 MELEMENT *mesh;
2096
2097 counter = 0;
2098
2099 for(i = 0; i < num_mesh; i++)
2100 for(mesh = Mlist[i]; mesh != NULL; mesh = mesh->mnext) {
2101 fprintf(fp, "%d\t%d\t%d\n", i+1, mesh->filindex+1, mesh->sign);
2102 counter++;
2103 }
2104
2105 if (counter != nnz)
2106 fprintf(stderr,"Internal Warning: nnz %d != counter %d\n",nnz,counter);
2107
2108 }
2109
dump_M_to_matlab(fp,Mlist,num_mesh,nnz,mname)2110 dump_M_to_matlab(fp, Mlist, num_mesh, nnz, mname)
2111 FILE *fp;
2112 MELEMENT **Mlist;
2113 int num_mesh, nnz;
2114 char *mname;
2115 {
2116
2117 double onerow[3];
2118 int i, counter;
2119 MELEMENT *mesh;
2120
2121 counter = 0;
2122
2123 for(i = 0; i < num_mesh; i++) {
2124 onerow[0] = i + 1;
2125 for(mesh = Mlist[i]; mesh != NULL; mesh = mesh->mnext) {
2126 onerow[1] = mesh->filindex + 1;
2127 onerow[2] = mesh->sign;
2128 savemat_mod(fp, machine+100, mname, nnz, 3, 0, onerow,
2129 (double *)NULL, counter++, 3);
2130 }
2131 }
2132
2133 if (counter != nnz)
2134 fprintf(stderr,"Internal Warning: nnz %d != counter %d\n",nnz,counter);
2135 }
2136
2137 /* this picks one node in each tree to be a ground node */
pick_ground_nodes(indsys)2138 pick_ground_nodes(indsys)
2139 SYS *indsys;
2140 {
2141 TREE *atree;
2142 SEGMENT *seg;
2143 seg_ptr segp;
2144 char type;
2145 NODES *node;
2146
2147 for(atree = indsys->trees; atree != NULL; atree = atree->next) {
2148 type = atree->loops->path->seg.type;
2149 if (type == NORMAL)
2150 node = getrealnode(((SEGMENT *)atree->loops->path->seg.segp)->node[0]);
2151 else
2152 node = getrealnode(((PSEUDO_SEG *)atree->loops->path->seg.segp)->node[0]);
2153
2154 if (node->index != -1) {
2155 fprintf(stderr,"huh? index == %d in pick_ground_node\n",node->index);
2156 exit(1);
2157 }
2158 node->index = 0;
2159 }
2160 }
2161
pick_subset(portlist,indsys)2162 int pick_subset(portlist, indsys)
2163 strlist *portlist;
2164 SYS *indsys;
2165 {
2166 strlist *oneport;
2167 EXTERNAL *ext;
2168 int counter;
2169
2170 if (portlist == NULL) {
2171 for(ext = indsys->externals; ext != NULL; ext=ext->next)
2172 ext->col_Yindex = ext->Yindex;
2173 return indsys->num_extern;
2174 }
2175
2176 for(ext = indsys->externals; ext != NULL; ext=ext->next)
2177 ext->col_Yindex = -1;
2178
2179 counter = 0;
2180 for(oneport = portlist; oneport != NULL; oneport = oneport->next) {
2181 ext = get_external_from_portname(oneport->str,indsys);
2182 if (ext == NULL) {
2183 fprintf(stderr,"Error: unknown portname %s\n",oneport->str);
2184 exit(1);
2185 }
2186
2187 ext->col_Yindex = counter;
2188 counter++;
2189 }
2190
2191 return counter;
2192 }
2193
2194 /* concatenates so that s1 = s1 + s2 + s3 + s4 */
concat4(s1,s2,s3,s4)2195 concat4(s1,s2,s3,s4)
2196 char *s1, *s2, *s3, *s4;
2197 {
2198 s1[0] = '\0';
2199 strcat(s1,s2);
2200 strcat(s1,s3);
2201 strcat(s1,s4);
2202 }
2203
2204