1 /* 2 Copyright (c) 1992 Massachusetts Institute of Technology, Cambridge, MA. 3 All rights reserved. 4 5 This Agreement gives you, the LICENSEE, certain rights and obligations. 6 By using the software, you indicate that you have read, understood, and 7 will comply with the terms. 8 9 Permission to use, copy and modify for internal, noncommercial purposes 10 is hereby granted. Any distribution of this program or any part thereof 11 is strictly prohibited without prior written consent of M.I.T. 12 13 Title to copyright to this software and to any associated documentation 14 shall at all times remain with M.I.T. and LICENSEE agrees to preserve 15 same. LICENSEE agrees not to make any copies except for LICENSEE'S 16 internal noncommercial use, or to use separately any portion of this 17 software without prior written consent of M.I.T. LICENSEE agrees to 18 place the appropriate copyright notice on any such copies. 19 20 Nothing in this Agreement shall be construed as conferring rights to use 21 in advertising, publicity or otherwise any trademark or the name of 22 "Massachusetts Institute of Technology" or "M.I.T." 23 24 M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. By 25 way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR 26 WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR 27 THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL 28 NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. 29 M.I.T. shall not be held liable for any liability nor for any direct, 30 indirect or consequential damages with respect to any claim by LICENSEE 31 or any third party on account of or arising from this Agreement or use 32 of this software. 33 */ 34 35 /* S. R. Whiteley (stevew@srware.com) adapted this code for 36 * superconductors 6/29/96. These modifications are active 37 * when the preprocessor variable SUPERCON is defined. 38 */ 39 40 /* # ***** sort to /src/header 41 # ***** */ 42 43 /* for NWS-3860 compatability */ 44 #ifdef NEWS 45 /* 46 ** Memory management functions (from stdlib.h (recursive incl => unusable)) 47 */ 48 extern char * calloc(); 49 extern char * malloc(); 50 extern char * realloc(); 51 #else 52 #include <stdlib.h> 53 #endif /* end if NEWS */ 54 #include <stdio.h> 55 #include <math.h> 56 #include <unistd.h> 57 58 /* fastcap data structures */ 59 #include "mulStruct.h" 60 61 /* execution time macros */ 62 #include "resusage.h" 63 64 /* time variables/structs */ 65 #ifndef _TIME_ /* if not on a Sun4 */ 66 #ifndef NEWS /* if not on a NWS-38XX */ 67 #include <time.h> 68 #endif 69 #endif 70 71 #define VERSION 3.0 72 73 /*********************************************************************** 74 macros for allocation with checks for NULL pntrs and 0 byte requests 75 - also keep an allocated memory count 76 - CALLOC() is used when the memory must be zeroed 77 its core should be either calloc() or ualloc() (not as fast but 78 more space efficient, no free list - uses sbrk() and never frees) 79 - MALLOC() used when memory can be anything 80 core should be malloc() or ualloc() 81 ***********************************************************************/ 82 /* #define CALCORE(NUM, TYPE) calloc((unsigned)(NUM),sizeof(TYPE)) */ 83 #define CALCORE(NUM, TYPE) ualloc((unsigned)(NUM)*sizeof(TYPE)) 84 /* #define MALCORE malloc */ 85 #define MALCORE ualloc 86 87 /* declare ualloc */ 88 extern char *ualloc(); 89 90 /* counts of memory usage by multipole matrix type */ 91 extern long memcount; 92 extern long memQ2M; 93 extern long memQ2L; 94 extern long memQ2P; 95 extern long memL2L; 96 extern long memM2M; 97 extern long memM2L; 98 extern long memM2P; 99 extern long memL2P; 100 extern long memQ2PD; 101 extern long memMSC; 102 extern long memIND; 103 104 #ifdef MATTDEBUG 105 extern long membins[1001]; 106 #endif 107 108 /* types of memory usage by multipole matrix type */ 109 #define AQ2M 0 110 #define AQ2L 1 111 #define AQ2P 2 112 #define AL2L 3 113 #define AM2M 4 114 #define AM2L 5 115 #define AM2P 6 116 #define AL2P 7 117 #define AQ2PD 8 118 #define AMSC 9 119 #define IND 10 120 121 #ifdef NO_SBRK 122 #define sbrk(x) 0L 123 #endif 124 125 #define DUMPALLOCSIZ \ 126 { \ 127 (void)fprintf(stderr, \ 128 "Total Memory Allocated: %ld kilobytes (brk = 0x%lx)\n", \ 129 memcount/1024, (long)sbrk(0)); \ 130 /* # ***** awked out for release */ \ 131 (void)fprintf(stderr, " Q2M matrix memory allocated: %7.ld kilobytes\n",\ 132 memQ2M/1024); \ 133 memcount = memQ2M; \ 134 (void)fprintf(stderr, " Q2L matrix memory allocated: %7.ld kilobytes\n",\ 135 memQ2L/1024); \ 136 memcount += memQ2L; \ 137 (void)fprintf(stderr, " Q2P matrix memory allocated: %7.ld kilobytes\n",\ 138 memQ2P/1024); \ 139 memcount += memQ2P; \ 140 (void)fprintf(stderr, " L2L matrix memory allocated: %7.ld kilobytes\n",\ 141 memL2L/1024); \ 142 memcount += memL2L; \ 143 (void)fprintf(stderr, " M2M matrix memory allocated: %7.ld kilobytes\n",\ 144 memM2M/1024); \ 145 memcount += memM2M; \ 146 (void)fprintf(stderr, " M2L matrix memory allocated: %7.ld kilobytes\n",\ 147 memM2L/1024); \ 148 memcount += memM2L; \ 149 (void)fprintf(stderr, " M2P matrix memory allocated: %7.ld kilobytes\n",\ 150 memM2P/1024); \ 151 memcount += memM2P; \ 152 (void)fprintf(stderr, " L2P matrix memory allocated: %7.ld kilobytes\n",\ 153 memL2P/1024); \ 154 memcount += memL2P; \ 155 (void)fprintf(stderr, " Q2PD matrix memory allocated: %7.ld kilobytes\n",\ 156 memQ2PD/1024); \ 157 memcount += memQ2PD; \ 158 (void)fprintf(stderr, " Miscellaneous mem. allocated: %7.ld kilobytes\n",\ 159 memMSC/1024); \ 160 memcount += memMSC; \ 161 (void)fprintf(stderr, " Inductance mem. allocated: %7.ld kilobytes\n",\ 162 memIND/1024); \ 163 memcount += memMSC; \ 164 (void)fprintf(stderr, " Total memory (check w/above): %7.ld kilobytes\n",\ 165 memcount/1024); \ 166 /* # ***** awked out for release */ \ 167 } 168 169 #define CALLOC(PNTR, NUM, TYPE, FLAG, MTYP) \ 170 { \ 171 if((NUM)*sizeof(TYPE)==0) \ 172 (void)fprintf(stderr, \ 173 "zero element request in file `%s' at line %d\n", \ 174 __FILE__, __LINE__); \ 175 else if(((PNTR)=(TYPE*)CALCORE(NUM, TYPE))==NULL) { \ 176 (void)fprintf(stderr, \ 177 "\nfastcap: out of memory in file `%s' at line %d\n", \ 178 __FILE__, __LINE__); \ 179 (void)fprintf(stderr, " (NULL pointer on %ld byte request)\n", \ 180 (NUM)*sizeof(TYPE)); \ 181 DUMPALLOCSIZ; \ 182 DUMPRSS; \ 183 (void)fflush(stderr); \ 184 (void)fflush(stdout); \ 185 if(FLAG == ON) exit(1); \ 186 } \ 187 else { \ 188 memcount += ((NUM)*sizeof(TYPE)); \ 189 if(MTYP == AQ2M) memQ2M += ((NUM)*sizeof(TYPE)); \ 190 else if(MTYP == AQ2L) memQ2L += ((NUM)*sizeof(TYPE)); \ 191 else if(MTYP == AQ2P) memQ2P += ((NUM)*sizeof(TYPE)); \ 192 else if(MTYP == AL2L) memL2L += ((NUM)*sizeof(TYPE)); \ 193 else if(MTYP == AM2M) memM2M += ((NUM)*sizeof(TYPE)); \ 194 else if(MTYP == AM2L) memM2L += ((NUM)*sizeof(TYPE)); \ 195 else if(MTYP == AM2P) memM2P += ((NUM)*sizeof(TYPE)); \ 196 else if(MTYP == AL2P) memL2P += ((NUM)*sizeof(TYPE)); \ 197 else if(MTYP == AQ2PD) memQ2PD += ((NUM)*sizeof(TYPE)); \ 198 else if(MTYP == AMSC) memMSC += ((NUM)*sizeof(TYPE)); \ 199 else if(MTYP == IND) memIND += ((NUM)*sizeof(TYPE)); \ 200 else { \ 201 (void)fprintf(stderr, "CALLOC: unknown memory type %d\n", MTYP); \ 202 exit(1); \ 203 } \ 204 /*if ((NUM)*sizeof(TYPE) >= 10000) \ 205 membins[1000] += 1; \ 206 else \ 207 membins[(NUM)*sizeof(TYPE)/10] += 1; */ \ 208 } \ 209 } 210 211 #define MALLOC(PNTR, NUM, TYPE, FLAG, MTYP) \ 212 { \ 213 if((NUM)*sizeof(TYPE)==0) \ 214 (void)fprintf(stderr, \ 215 "zero element request in file `%s' at line %d\n", \ 216 __FILE__, __LINE__); \ 217 else if(((PNTR)=(TYPE*)MALCORE((unsigned)((NUM)*sizeof(TYPE))))==NULL) { \ 218 (void)fprintf(stderr, \ 219 "\nfastcap: out of memory in file `%s' at line %d\n", \ 220 __FILE__, __LINE__); \ 221 (void)fprintf(stderr, " (NULL pointer on %ld byte request)\n", \ 222 (NUM)*sizeof(TYPE)); \ 223 DUMPALLOCSIZ; \ 224 DUMPRSS; \ 225 (void)fflush(stderr); \ 226 (void)fflush(stdout); \ 227 if(FLAG == ON) exit(1); \ 228 } \ 229 else { \ 230 memcount += ((NUM)*sizeof(TYPE)); \ 231 if(MTYP == AQ2M) memQ2M += ((NUM)*sizeof(TYPE)); \ 232 else if(MTYP == AQ2L) memQ2L += ((NUM)*sizeof(TYPE)); \ 233 else if(MTYP == AQ2P) memQ2P += ((NUM)*sizeof(TYPE)); \ 234 else if(MTYP == AL2L) memL2L += ((NUM)*sizeof(TYPE)); \ 235 else if(MTYP == AM2M) memM2M += ((NUM)*sizeof(TYPE)); \ 236 else if(MTYP == AM2L) memM2L += ((NUM)*sizeof(TYPE)); \ 237 else if(MTYP == AM2P) memM2P += ((NUM)*sizeof(TYPE)); \ 238 else if(MTYP == AL2P) memL2P += ((NUM)*sizeof(TYPE)); \ 239 else if(MTYP == AQ2PD) memQ2PD += ((NUM)*sizeof(TYPE)); \ 240 else if(MTYP == AMSC) memMSC += ((NUM)*sizeof(TYPE)); \ 241 else if(MTYP == IND) memIND += ((NUM)*sizeof(TYPE)); \ 242 else { \ 243 (void)fprintf(stderr, "MALLOC: unknown memory type %d\n", MTYP); \ 244 exit(1); \ 245 } \ 246 /*if ((NUM)*sizeof(TYPE) >= 10000) \ 247 membins[1000] += 1; \ 248 else \ 249 membins[(NUM)*sizeof(TYPE)/10] += 1; */ \ 250 } \ 251 } 252 253 /***************************************************************************** 254 255 misc. global macros 256 257 *****************************************************************************/ 258 #define NOT ! 259 #define ABORT() \ 260 { (void)fflush(stdout); \ 261 (void)fprintf(stderr, "FastCap: panic in file `%s' at line %d.\n",\ 262 __FILE__, __LINE__); \ 263 (void)fflush(stderr); \ 264 abort(); \ 265 } 266 267 #define ASSERT(condition) if(NOT(condition)) ABORT() 268 269 #define INNER(pap,p,ap,size) for(pap=0.0,i=1; i<=size; i++) pap += p[i]*ap[i]; 270 271 #ifndef MAX 272 #define MAX(A,B) ( (A) > (B) ? (A) : (B) ) 273 #endif 274 275 #ifndef MIN 276 #define MIN(A,B) ( (A) > (B) ? (B) : (A) ) 277 #endif 278 279 #define ABS(A) ( ( (A) > 0 ) ? (A) : (-(A)) ) 280 281 #define VCOPY(A, B) A[0] = B[0]; A[1] = B[1]; A[2] = B[2]; 282 283 #define TRUE 1 284 #define FALSE 0 285 286 #define ON 1 287 #define OFF 0 288 289 #define LAST 2 290 #define ALL 2 291 292 #ifndef M_PI 293 /* pi constant included here since won't be in ANSI C */ 294 #define M_PI 3.1415926535897931160E0 /*Hex 2^ 1 * 1.921FB54442D18 */ 295 #endif 296 297 #define E_0 8.854187818E-12 /* epsilon0 +- .000000071E-12 F/m */ 298 #define FPIEPS 4.0*M_PI*E_0 /* 4 pi times the dielectric permittivity, 299 free-space permittivity is the default, 300 units are F/m - all dimensions in meters */ 301 302 /* flags in chkList() in mulDisplay.c (chks direct, local or eval cube lsts) */ 303 #define DIRECT 0 304 #define LOCAL 1 305 #define EVAL 3 306 307 /* types of surfaces */ 308 #define CONDTR 0 /* conductor surfaces */ 309 #define DIELEC 1 /* dielectric interface surface */ 310 #define BOTH 3 /* dielectric i/f w/very thin cond on it */ 311 312 /* used in input routines */ 313 #define MAXCON 10000 /* assumes never more conductors than this */ 314 315 /* used in ps file dump */ 316 #define OPEN 0 /* open ps file, print hdr, ignore row/col */ 317 #define CLOSE 1 /* print trailer, close ps file */ 318 #define UPDATE 2 /* => add 2 dots for this row and col */ 319 320 /* divided difference distances, see electric.c */ 321 #define HPOS (1e-6*cur_panel->max_diag) /* h in positive normal dir */ 322 #define HNEG HPOS /* h for divided difference in neg nrml dir */ 323 324 /* level set mode, see placeq, mulSetup.c and input.c */ 325 #define ONELES 2 /* => auto set levs to 1 up fr fully exact */ 326 327 /* expansion moment index generating macros (see mulMulti.c, mulLocal.c) */ 328 #define CINDEX(N, M) ( (M) + ((N)*((N)+1))/2 ) 329 #define SINDEX(N, M, CTERMS) ( (CTERMS) + (M) + ((N)*((N)+1))/2 - ((N)+1) ) 330 331 /* used in get_kill_num_list and routines it calls */ 332 #define NOTUNI -1 333 #define NOTFND -2 334 335 /*********************************************************************** 336 337 configuration and debug flags 338 339 ***********************************************************************/ 340 341 /* types of downward/eval passes */ 342 #define NOLOCL 0 /* multipoles evaluated directly - no locals */ 343 #define NOSHFT 1 /* multis to locals w/o local2local shifts */ 344 #define GRENGD 3 /* full Greengard downward pass/eval */ 345 346 /* types of iterative methods (values of ITRTYP below) */ 347 #define GCR 0 /* GCR with single (not block) vector iters */ 348 #define GMRES 1 /* GMRES with vector iterates */ 349 350 /* types of finite elements (NOTE: only const. chg den. panels implemented) */ 351 #define CONST 0 /* constant charge density on panels */ 352 #define AFFINE 1 353 #define QUADRA 2 354 355 /* types of weighted residuals methods (NOTE: only collocation implemented) */ 356 #define COLLOC 0 /* point collocation */ 357 #define SUBDOM 1 /* subdomain collocation */ 358 #define GALKIN 2 /* Galerkin */ 359 360 /* types of preconditioners. */ 361 #define NONE 0 362 #define BD 1 /* Block diagonal (not set up for dielecs). */ 363 #define OL 2 /* OverLap */ 364 365 /* Discretization Configuration */ 366 #define WRMETH COLLOC /* weighted res meth type (COLLOC only now) */ 367 #define ELTYPE CONST /* finite element type (CONST only now) */ 368 /* Multipole Configuration */ 369 #define DNTYPE GRENGD /* type of downward/eval pass - see above */ 370 #define MULTI ON /* ON=> add in multipole contribution to P*q */ 371 #define RADINTER ON /* ON=> Parent level multis in interlist. */ 372 #define NNBRS 2 /* Distance to consider a nearest nbr. */ 373 #define ADAPT ON /* ON=> use adaptive algorithm */ 374 #define OPCNT OFF /* Counts the Matrix-Vector multiply ops. */ 375 #define DEFORD 2 /* default expansion order */ 376 #define MAXORDER 6 /* Maximum expansion order (sets ary sizes) */ 377 #define MAXDEP 20 /* maximum partitioning depth */ 378 #define NUMDPT 2 /* num pnts for ea dielec panel (2 or 3) */ 379 #define SKIPQD OFF /* ON => skip dielec panel chg in E eval */ 380 /* Linear System Solution Configuration */ 381 #define ITRTYP GMRES /* type of iterative method */ 382 #define PRECOND NONE /* NONE=> no preconditioner OL=> use prec. */ 383 #define DIRSOL OFF /* ON=> solve Pq=psi by Gaussian elim. */ 384 #define EXPGCR OFF /* ON=> do explicit full P*q products */ 385 #define ABSTOL 0.01 /* iterations until ||res||inf < ABSTOL */ 386 #define MAXITER size /* max num iterations ('size' => # panels) */ 387 #define EXRTSH 0.9 /* exact/ttl>EXRTSH for lev => make last lev */ 388 #define SUPERCON ON /* handle superconductive conductors */ 389 /* (add any new configuration flags to dumpConfig() in mulDisplay.c) */ 390 391 /* Output Format Configuration */ 392 #define MKSDAT ON /* ON=> dump symmetrized, MKS units cap mat */ 393 #define CMDDAT ON /* ON=> dump command line info to output */ 394 #define RAWDAT OFF /* ON=> dump unsymm, Gaussian units cap mat */ 395 #define ITRDAT OFF /* ON=> dump residuals for every iteration */ 396 #define TIMDAT ON /* ON=> dump time and memory usage numbers */ 397 #define CFGDAT OFF /* ON=> dump configuration flags to output */ 398 #define MULDAT OFF /* ON=> dump brief multipole setup info */ 399 #define DISSYN OFF /* ON=> display synopsis of cubes in lists */ 400 #define DMTCNT OFF /* ON=> display xform matrix counts by level */ 401 #define DISSRF OFF /* ON=> display input surface information */ 402 #define NAMDAT OFF /* ON=> dump conductor names */ 403 #define CAPVEW ON /* ON=> enable ps file dumps of geometry */ 404 405 /* display of transformation matrices */ 406 #define DISQ2M OFF /* ON=> display Q2M matrices when built */ 407 #define DISM2M OFF /* ON=> display M2M matrices when built */ 408 #define DISM2P OFF /* ON=> display M2P matrices when built */ 409 #define DISL2P OFF /* ON=> display L2P matrices when built */ 410 #define DISQ2P OFF /* ON=> display Q2P matrices when built */ 411 #define DSQ2PD OFF /* ON=> display Q2PDiag matrices > build */ 412 #define DISQ2L OFF /* ON=> display Q2L matrices when built */ 413 #define DISM2L OFF /* ON=> display M2L matrices when built */ 414 #define DISL2L OFF /* ON=> display L2L matrices when built */ 415 #define DALQ2M OFF /* ON=> display all Q2M matrix build steps */ 416 #define DALM2P OFF /* ON=> display all M2P matrix build steps */ 417 #define DALL2P OFF /* ON=> display all L2P matrix build steps */ 418 #define DALQ2L OFF /* ON=> display all Q2L matrix build steps */ 419 /* display of other intermediate results */ 420 #define DUPVEC OFF /* ON=> display lev 1 upward pass vectors */ 421 #define DISFAC OFF /* ON=> display factorial fractions in M2L */ 422 #define DPSYSD OFF /* ON=> display system after direct build */ 423 #define DILIST OFF /* ON=> display interaction lists */ 424 #define DMPELE OFF /* ON=> display electric flux densities */ 425 #define DMPCHG OFF /* ON=> display all charge vector iterates 426 LAST=> display final charge vector */ 427 /* misc debug */ 428 #define CKDLST OFF /* ON=> check direct list, prnt msg if bad */ 429 #define DMPREC OFF /* ON=> dump P and Ctil to matlab file */ 430 #define CKCLST OFF /* ON=> check charge list, prnt msg if bad */ 431 #define DUMPPS OFF /* ON=> dump ps file w/mulMatDirect calcp's 432 ALL=> dump adaptive alg calcp's as well */ 433 #define DPCOMP OFF /* ON=> dump prec pts before&aft compression */ 434 #define DPDDIF OFF /* ON=> dump divided difference components */ 435 #define CHKDUM OFF /* ON=> print msg if dummy list inconsistent */ 436 #define JACDBG OFF /* ON=> print random Jacob debug messages */ 437 /* blkDirect.c related flags - used only when DIRSOL == ON || EXPGCR == ON */ 438 #define MAXSIZ 0 /* any more tiles than this uses matrix on disk 439 for DIRSOL == ON or EXPGCR == ON */ 440