This is gck_buildsys.c in view mode; [Download] [Up]
/* ******************************************************************** */
/* * GCK * */
/* * A Circuit Simulation Program * */
/* * by Tanju Cataltepe * */
/* * (c) Copyright 1989 * */
/* ******************************************************************** */
/* (c) Copyright 1989, Tanju Cataltepe */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#ifdef Macintosh
#include <console.h>
#include <Files.h>
#else
#include <sys/types.h>
#include <sys/stat.h>
#endif
#include <math.h>
#include <float.h>
#include "gck.h"
#include "gck_vars.h"
#include "gck_protos.h"
#include "gck_errors.h"
void gaussj(nv)
long nv;
{
/* RdH Debug */
long nsn; /* node search number */
/*Solves the equation Mx = Ny + Bs for x.
M -> ?
N -> (M_inv)N
B -> (M_inv)B
*/
double *indxc;
double *indxr;
double *ipiv;
double big, dum, pivinv;
long i, icol, irow, j, k, l, ll;
indxc = (double *) c_malloc(nv * sizeof(double));
indxr = (double *) c_malloc(nv * sizeof(double));
ipiv = (double *) c_malloc(nv * sizeof(double));
/* Starting Gaussj M */
for (j = 0; j < nv; j++) {
indxc[j] = 0.0; /* new RdH */
indxr[j] = 0.0; /* new RdH */
ipiv[j] = 0.0;
}
for (i = 0; i < nv; i++) {
big = 0.0;
for (j = 1; j <= nv; j++) {
if (ipiv[j - 1] != 1) {
for (k = 1; k <= nv; k++) {
if (ipiv[k - 1] == 0) {
if (fabs(M[j - 1][k - 1]) >= big) {
big = fabs(M[j - 1][k - 1]);
irow = j;
icol = k;
}
} else if (ipiv[k - 1] > 1)
errorReport(ERROR_GAUSSIAN, progName, ERR_VOID, 0);
}
}
}
ipiv[icol - 1]++;
if (irow != icol) {
for (l = 0; l < nv; l++) {
dum = M[irow - 1][l];
M[irow - 1][l] = M[icol - 1][l];
M[icol - 1][l] = dum;
dum = N[irow - 1][l];
N[irow - 1][l] = N[icol - 1][l];
N[icol - 1][l] = dum;
dum = B[irow - 1][l];
B[irow - 1][l] = B[icol - 1][l];
B[icol - 1][l] = dum;
}
}
indxr[i] = irow;
indxc[i] = icol;
if (M[icol - 1][icol - 1] == 0) {
for (nsn = 0; nsn < ID_nindex_size; nsn++) {
if (ID_nindex[nsn] == icol)
errorReport(ERROR_GAUSSIAN, progName, ERR_NODE, nsn);
}
errorReport(ERROR_GAUSSIAN, progName, ERR_NODE, 0);
}
pivinv = 1 / M[icol - 1][icol - 1];
M[icol - 1][icol - 1] = 1.0;
for (l = 0; l < nv; l++) {
M[icol - 1][l] *= pivinv;
N[icol - 1][l] *= pivinv;
B[icol - 1][l] *= pivinv;
}
for (ll = 0; ll < nv; ll++) {
if (ll + 1 != icol) {
dum = M[ll][icol - 1];
M[ll][icol - 1] = 0.0;
for (l = 0; l < nv; l++) {
M[ll][l] -= M[icol - 1][l] * dum;
N[ll][l] -= N[icol - 1][l] * dum;
B[ll][l] -= B[icol - 1][l] * dum;
}
}
}
}
/* End Gaussj */
free(indxc);
free(indxr);
free(ipiv);
}
void BuildSystem()
{
long i, j, k;
long matrix_loop;
systemRec *tmpSystem;
branchRec *currentBranch, *controlBranch;
long bn, n1, n2, nc1, nc2;
double val;
long FORLIM, FORLIM1, FORLIM2;
FILE *fp;
aWord compiled_data_filename;
boolean writing_to_disk;
writing_to_disk = TRUE;
if (cmdline_symbols == TRUE)
writing_to_disk = FALSE;
strcpy(compiled_data_filename,rfile.filename);
strcat(compiled_data_filename,".gck");
if ((cmdline_symbols == TRUE) &&
(flag_compiler == COMPILER_COMPILE_AND_QUIT)) {
remove(compiled_data_filename);
errorReport(ERROR_SYMBOLS_ON_CQ,0,ERR_VOID,0);
}
if (writing_to_disk) {
fp = fopen(compiled_data_filename,"wb");
if (fp == NULL)
errorReport(ERROR_COMPILER_FILE_COULDNT_OPEN,compiled_data_filename,ERR_VOID,0);
fwrite(&NoOfTicks,sizeof(long),1,fp);
fwrite(&NoOfVars,sizeof(long),1,fp);
}
#ifdef Macintosh
if (PRINT_HIGHLIGHTS)
fprintf(stderr,"*Mac-Specific Warning : The following Procedure may take a while!\n");
#endif
if (PRINT_PAINFUL)
fprintf(stderr,"Allocating Arrays\n");
M = (double **) c_malloc(NoOfVars * sizeof (double *));
N = (double **) c_malloc(NoOfVars * sizeof (double *));
B = (double **) c_malloc(NoOfVars * sizeof (double *));
#ifdef Macintosh
if (PRINT_HIGHLIGHTS)
fprintf(stderr,"Allocating %ld for M\n",NoOfVars);
#endif
for (i = 0; i < NoOfVars; i++) {
#ifdef Macintosh
if (PRINT_PAINFUL)
fprintf(stderr,"M %5ld\n",i);
#endif
M[i] = (double *) c_malloc(NoOfVars * sizeof (double));
}
#ifdef Macintosh
if (PRINT_HIGHLIGHTS)
fprintf(stderr,"Allocating %ld for N\n",NoOfVars);
#endif
for (i = 0; i < NoOfVars; i++) {
#ifdef Macintosh
if (PRINT_PAINFUL)
fprintf(stderr,"N %5ld\n",i);
#endif
N[i] = (double *) c_malloc(NoOfVars * sizeof (double));
}
#ifdef Macintosh
if (PRINT_HIGHLIGHTS)
fprintf(stderr,"Allocating %ld for B\n",NoOfVars);
#endif
for (i = 0; i < NoOfVars; i++) {
#ifdef Macintosh
if (PRINT_PAINFUL)
fprintf(stderr,"B %5ld\n",i);
#endif
B[i] = (double *) c_malloc(NoOfVars * sizeof (double));
}
#ifdef Macintosh
if (PRINT_HIGHLIGHTS)
fprintf(stderr,"Complete.\n");
#endif
FORLIM = NoOfTicks;
for (i = 1; i <= FORLIM; i++) {
GetFromSystem(&tmpSystem, i);
FORLIM1 = NoOfVars;
for (k = 0; k < FORLIM1; k++) {
FORLIM2 = NoOfVars;
for (j = 0; j < FORLIM2; j++) {
M[k][j] = 0.0;
N[k][j] = 0.0;
B[k][j] = 0.0;
}
}
currentBranch = baseOfBranch;
while (currentBranch != NULL) {
if (currentBranch->element != subckt) {
/*Connectivity description*/
/* RdH */
/* bn : branch number */
/* n1 : node 1 (from) */
/* n2 : node 2 (to) */
bn = GetID(currentBranch->ID, 1L);
n1 = GetID(currentBranch->ID, 2L);
n2 = GetID(currentBranch->ID, 3L);
if (n1 != 0)
M[n1 - 1][bn - 1] = 1.0;
if (n2 != 0)
M[n2 - 1][bn - 1] = -1.0;
/*branch characteristics*/
if (((1L << ((long)currentBranch->element)) & fourTerminal) != 0) {
nc1 = GetID(currentBranch->ID, 4L);
nc2 = GetID(currentBranch->ID, 5L);
}
if (((1L << ((long)currentBranch->element)) &
((1L << ((long)ccvs)) | (1L << ((long)cccs)))) != 0) {
GetBranchByName(&controlBranch, currentBranch->cname);
nc1 = GetID(controlBranch->ID, 1L);
}
if (((1L << ((long)currentBranch->element)) & ((1L << ((long)quantizer)) |
(1L << ((long)vs)) | (1L << ((long)cs)))) != 0)
val = 1.0;
else if (((1L << ((long)currentBranch->element)) &
(1L << ((long)switch_))) != 0)
val = GetValue(currentBranch->value, i);
else if (((1L << ((long)currentBranch->element)) & digitals) == 0 &&
((1L << ((long)currentBranch->element)) & logicals) == 0)
val = GetValue(currentBranch->value, 1L);
switch (currentBranch->element) {
case resistor:
M[bn - 1][bn - 1] = -val;
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
break;
case capacitor:
val /= deltaT;
M[bn - 1][bn - 1] = -1.0;
if (n1 != 0) {
M[bn - 1][n1 - 1] = val;
N[bn - 1][n1 - 1] = val;
}
if (n2 != 0) {
M[bn - 1][n2 - 1] = -val;
N[bn - 1][n2 - 1] = -val;
}
break;
case inductor:
val /= deltaT;
M[bn - 1][bn - 1] = -val;
N[bn - 1][bn - 1] = -val;
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
break;
case vcvs:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
if (nc1 != 0)
M[bn - 1][nc1 - 1] = -val;
if (nc2 != 0)
M[bn - 1][nc2 - 1] = val;
break;
case vccs:
M[bn - 1][bn - 1] = 1.0;
if (nc1 != 0)
M[bn - 1][nc1 - 1] = -val;
if (nc2 != 0)
M[bn - 1][nc2 - 1] = val;
break;
case ccvs:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
M[bn - 1][nc1 - 1] = -val;
break;
case cccs:
M[bn - 1][bn - 1] = 1.0;
M[bn - 1][nc1 - 1] = -val;
break;
case switch_:
if (val == 0)
M[bn - 1][bn - 1] = 1.0;
else {
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
}
break;
case delay:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case adder:
if (n1 != 0)
M[bn - 1][n1 - 1] = -1.0;
if (nc1 != 0)
M[bn - 1][nc1 - 1] = GetValue(currentBranch->value, 1L);
if (nc2 != 0)
M[bn - 1][nc2 - 1] = GetValue(currentBranch->value, 2L);
break;
case multiplier:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case opamp:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
if (nc1 != 0)
M[bn - 1][nc1 - 1] = -val;
if (nc2 != 0)
M[bn - 1][nc2 - 1] = val;
break;
case quantizer:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case vs:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case cs:
M[bn - 1][bn - 1] = 1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case l_general:
if (n1 != 0)
M[bn - 1][n1 - 1] = 1.0;
if (n2 != 0)
M[bn - 1][n2 - 1] = -1.0;
B[bn - 1][bn - 1] = 1.0;
break;
case subckt:
/* blank case */
break;
default:
errorReport(ERROR_GAUSSIAN_UNKNOWN_ELEMENT, currentBranch->name, ERR_VOID, 0);
break;
}
}
currentBranch = currentBranch->next;
}
/*Do Gaussian elimination on Mx=Ny+Bs simultaneously*/
if (runDiag)
Diagnostic(); /* No longer pass parameters...arrays are global */
gaussj(NoOfVars);
/* Initialize new tmpSystem records */
FORLIM1 = NoOfVars;
for (k = 0; k < FORLIM1; k++) {
FORLIM2 = NoOfVars;
for (j = 0; j < FORLIM2; j++) {
tmpSystem->N[k][j] = N[k][j];
tmpSystem->B[k][j] = B[k][j];
}
if (writing_to_disk) {
fwrite((tmpSystem->N[k]),sizeof(double),(NoOfVars),fp);
fwrite((tmpSystem->B[k]),sizeof(double),(NoOfVars),fp);
}
}
if (PRINT_PAINFUL)
fprintf(stderr,"Finished Writing Circuit Phase %8ld to Disk.\n",i);
}
/* free memory for M, N, B */
if (writing_to_disk)
fclose(fp);
if (cmdline_symbols == TRUE)
remove(compiled_data_filename);
for (i = 0; i < NoOfVars; i++) {
free(M[i]);
free(N[i]);
free(B[i]);
}
free(M);
free(N);
free(B);
}
void BuildQuickSystem()
{
long i, j, k;
long NqTicks, NqVars;
systemRec *tmpSystem;
FILE *fp;
aWord compiled_data_filename;
strcpy(compiled_data_filename,rfile.filename);
strcat(compiled_data_filename,".gck");
fp = fopen(compiled_data_filename,"rb");
if (fp == NULL)
errorReport(ERROR_COMPILER_FILE_NOT_FOUND,compiled_data_filename,ERR_VOID,0);
fread(&NqTicks,sizeof(long),1,fp);
fread(&NqVars,sizeof(long),1,fp);
if (NqTicks != NoOfTicks)
errorReport(ERROR_COMPILER_TICKS_INCONSISTENT,0,ERR_VOID,0);
if (NqVars != NoOfVars)
errorReport(ERROR_COMPILER_VARS_INCONSISTENT,0,ERR_VOID,0);
for (i = 1; i <= NoOfTicks; i++) {
GetFromSystem(&tmpSystem, i);
for (k = 0; k < NoOfVars; k++) {
fread((tmpSystem->N[k]),sizeof(double),(NoOfVars),fp);
fread((tmpSystem->B[k]),sizeof(double),(NoOfVars),fp);
}
}
fclose(fp);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.