diff --git a/src/spicelib/devices/cpl/cplsetup.c b/src/spicelib/devices/cpl/cplsetup.c index 76865084e..786a35591 100644 --- a/src/spicelib/devices/cpl/cplsetup.c +++ b/src/spicelib/devices/cpl/cplsetup.c @@ -57,7 +57,7 @@ static double C_m[MAX_DIM][MAX_DIM]; static double length; static double TAU[MAX_DIM]; -static double A[MAX_DIM][2*MAX_DIM]; +static double A[MAX_DIM][2 * MAX_DIM]; static double frequency[MAX_DEG]; @@ -65,11 +65,11 @@ static double Si[MAX_DIM][MAX_DIM]; static double Si_1[MAX_DIM][MAX_DIM]; /* MacLaurin Series */ -static double *SiSv_1[MAX_DIM][MAX_DIM]; -static double *Sip[MAX_DIM][MAX_DIM]; -static double *Si_1p[MAX_DIM][MAX_DIM]; -static double *Sv_1p[MAX_DIM][MAX_DIM]; -static double *W[MAX_DIM]; +static double* SiSv_1[MAX_DIM][MAX_DIM]; +static double* Sip[MAX_DIM][MAX_DIM]; +static double* Si_1p[MAX_DIM][MAX_DIM]; +static double* Sv_1p[MAX_DIM][MAX_DIM]; +static double* W[MAX_DIM]; static Mult_Out IWI[MAX_DIM][MAX_DIM]; static Mult_Out IWV[MAX_DIM][MAX_DIM]; @@ -80,7 +80,7 @@ static double Scaling_F2; /* misc.c match */ static void new_memory(int, int, int); -static double *vector(int, int); +static double* vector(int, int); static void free_vector(double*, int, int); static void polint(double*, double*, int, double, double*, double*); static int match(int, double*, double*, double*); @@ -88,7 +88,7 @@ static int match(int, double*, double*, double*); static int Gaussian_Elimination2(int, int); static void eval_Si_Si_1(int, double); static void loop_ZY(int, double); -static void poly_matrix(double *A[MAX_DIM][MAX_DIM], int dim, int deg); +static void poly_matrix(double* A[MAX_DIM][MAX_DIM], int dim, int deg); /* static int checkW(double*, double); */ static void poly_W(int, int); static void eval_frequency(int, int); @@ -102,11 +102,11 @@ static int ReadCpL(CPLinstance*, CKTcircuit*); /* mult */ static void mult_p(double*, double*, double*, int, int, int); -static void matrix_p_mult(double *A[MAX_DIM][MAX_DIM], - double *D[MAX_DIM], - double *B[MAX_DIM][MAX_DIM], - int dim, int deg, int deg_o, - Mult_Out X[MAX_DIM][MAX_DIM]); +static void matrix_p_mult(double* A[MAX_DIM][MAX_DIM], + double* D[MAX_DIM], + double* B[MAX_DIM][MAX_DIM], + int dim, int deg, int deg_o, + Mult_Out X[MAX_DIM][MAX_DIM]); static double approx_mode(double*, double*, double); static double eval2(double, double, double, double); static int get_c(double, double, double, double, double, double, double, double*, double*); @@ -120,7 +120,7 @@ static NODE* insert_node(char*); static NDnamePt insert_ND(char*, NDnamePt*); static NODE* NEW_node(void); static NDnamePt ndn_btree; -static NODE *node_tab; +static NODE* node_tab; #define epsi_mult 1e-28 /* diag */ @@ -132,468 +132,469 @@ static void diag(int); static int rotate(int, int, int); #define epsi 1.0e-16 -static char *message = "tau of coupled lines is larger than max time step"; +static char* message = "tau of coupled lines is larger than max time step"; /* ARGSUSED */ int -CPLsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *state) +CPLsetup(SMPmatrix* matrix, GENmodel* inModel, CKTcircuit* ckt, int* state) { - CPLmodel *model = (CPLmodel *)inModel; - CPLinstance *here; - CKTnode *tmp, *node; - int error, m, p; - char **branchname; - int noL; - - NG_IGNORE(state); - - /* hash table for local gc */ - mem_init(); - - /* loop through all the models */ - for( ; model != NULL; model = CPLnextModel(model)) { - - if (!model->Rmgiven) { - SPfrontEnd->IFerrorf (ERR_FATAL, - "model %s: lossy line series resistance not given", model->CPLmodName); - return(E_BADPARM); - } - if (!model->Gmgiven) { - SPfrontEnd->IFerrorf (ERR_FATAL, - "model %s: lossy line parallel conductance not given", model->CPLmodName); - return(E_BADPARM); - } - if (!model->Lmgiven) { - SPfrontEnd->IFerrorf (ERR_FATAL, - "model %s: lossy line series inductance not given", model->CPLmodName); - return (E_BADPARM); - } - if (!model->Cmgiven) { - SPfrontEnd->IFerrorf (ERR_FATAL, - "model %s: lossy line parallel capacitance not given", model->CPLmodName); - return (E_BADPARM); - } - if (!model->lengthgiven) { - SPfrontEnd->IFerrorf (ERR_FATAL, - "model %s: lossy line length must be given", model->CPLmodName); - return (E_BADPARM); - } - - /* loop through all the instances of the model */ - for (here = CPLinstances(model); here != NULL ; - here=CPLnextInstance(here)) { - - if (!here->CPLlengthGiven) - here->CPLlength=0.0; - - /* macro to make elements with built in test for out of memory */ + CPLmodel* model = (CPLmodel*)inModel; + CPLinstance* here; + CKTnode* tmp, * node; + int error, m, p; + char** branchname; + int noL; + + NG_IGNORE(state); + + /* hash table for local gc */ + mem_init(); + + /* loop through all the models */ + for (; model != NULL; model = CPLnextModel(model)) { + + if (!model->Rmgiven) { + SPfrontEnd->IFerrorf(ERR_FATAL, + "model %s: lossy line series resistance not given", model->CPLmodName); + return(E_BADPARM); + } + if (!model->Gmgiven) { + SPfrontEnd->IFerrorf(ERR_FATAL, + "model %s: lossy line parallel conductance not given", model->CPLmodName); + return(E_BADPARM); + } + if (!model->Lmgiven) { + SPfrontEnd->IFerrorf(ERR_FATAL, + "model %s: lossy line series inductance not given", model->CPLmodName); + return (E_BADPARM); + } + if (!model->Cmgiven) { + SPfrontEnd->IFerrorf(ERR_FATAL, + "model %s: lossy line parallel capacitance not given", model->CPLmodName); + return (E_BADPARM); + } + if (!model->lengthgiven) { + SPfrontEnd->IFerrorf(ERR_FATAL, + "model %s: lossy line length must be given", model->CPLmodName); + return (E_BADPARM); + } + + /* loop through all the instances of the model */ + for (here = CPLinstances(model); here != NULL; + here = CPLnextInstance(here)) { + + if (!here->CPLlengthGiven) + here->CPLlength = 0.0; + + /* macro to make elements with built in test for out of memory */ #define TSTALLOC(ptr,first,second) \ do { if((here->ptr = SMPmakeElt(matrix, here->first, here->second)) == NULL){\ return(E_NOMEM);\ } } while(0) - noL = here->dimension; - - here->CPLposNodes = TMALLOC(int, noL); - here->CPLnegNodes = TMALLOC(int, noL); - here->CPLibr1 = TMALLOC(int, noL); - here->CPLibr2 = TMALLOC(int, noL); - - VECTOR_ALLOC(double, here->CPLibr1Ibr1Ptr, noL); - VECTOR_ALLOC(double, here->CPLibr2Ibr2Ptr, noL); - VECTOR_ALLOC(double, here->CPLposIbr1Ptr, noL); - VECTOR_ALLOC(double, here->CPLnegIbr2Ptr, noL); - VECTOR_ALLOC(double, here->CPLposPosPtr, noL); - VECTOR_ALLOC(double, here->CPLnegNegPtr, noL); - VECTOR_ALLOC(double, here->CPLnegPosPtr, noL); - VECTOR_ALLOC(double, here->CPLposNegPtr, noL); - - MATRIX_ALLOC(double, here->CPLibr1PosPtr, noL, noL); - MATRIX_ALLOC(double, here->CPLibr2NegPtr, noL, noL); - MATRIX_ALLOC(double, here->CPLibr1NegPtr, noL, noL); - MATRIX_ALLOC(double, here->CPLibr2PosPtr, noL, noL); - MATRIX_ALLOC(double, here->CPLibr1Ibr2Ptr, noL, noL); - MATRIX_ALLOC(double, here->CPLibr2Ibr1Ptr, noL, noL); - - - branchname = TMALLOC(char *, here->dimension); - if (! here->CPLibr1Given) { - for (m = 0; m < here->dimension; m++) { - branchname[m] = TMALLOC(char, MAX_STRING); - sprintf(branchname[m], "branch1_%d", m); - error = - CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]); - if (error) return (error); - here->CPLibr1[m] = tmp->number; - tfree(branchname[m]); - } - here->CPLibr1Given = 1; - } - tfree(branchname); - branchname = TMALLOC(char *, here->dimension); - if (! here->CPLibr2Given) { - for (m = 0; m < here->dimension; m++) { - branchname[m] = TMALLOC(char, MAX_STRING); - sprintf(branchname[m], "branch2_%d", m); - error = - CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]); - if (error) return (error); - here->CPLibr2[m] = tmp->number; - tfree(branchname[m]); - } - here->CPLibr2Given = 1; - } - tfree(branchname); - - for (m = 0; m < here->dimension; m++) { - for (node = ckt->CKTnodes; node; node = node->next) { - if (strcmp(here->in_node_names[m], - node->name) == 0) { - here->CPLposNodes[m] = node->number; - } - } - } - for (m = 0; m < here->dimension; m++) { - for (node = ckt->CKTnodes; node; node = node->next) { - if (strcmp(here->out_node_names[m], - node->name) == 0) { - here->CPLnegNodes[m] = node->number; - } - } - } - - for (m = 0; m < here->dimension; m++) { - TSTALLOC(CPLibr1Ibr1Ptr[m],CPLibr1[m],CPLibr1[m]); - TSTALLOC(CPLibr2Ibr2Ptr[m],CPLibr2[m],CPLibr2[m]); - TSTALLOC(CPLposIbr1Ptr[m],CPLposNodes[m],CPLibr1[m]); - TSTALLOC(CPLnegIbr2Ptr[m],CPLnegNodes[m],CPLibr2[m]); - TSTALLOC(CPLposPosPtr[m],CPLposNodes[m],CPLposNodes[m]); - TSTALLOC(CPLnegNegPtr[m],CPLnegNodes[m],CPLnegNodes[m]); - TSTALLOC(CPLnegPosPtr[m],CPLnegNodes[m],CPLposNodes[m]); - TSTALLOC(CPLposNegPtr[m],CPLposNodes[m],CPLnegNodes[m]); - - for (p = 0; p < here->dimension; p++) { - - TSTALLOC(CPLibr1PosPtr[m][p],CPLibr1[m],CPLposNodes[p]); - TSTALLOC(CPLibr2NegPtr[m][p],CPLibr2[m],CPLnegNodes[p]); - TSTALLOC(CPLibr1NegPtr[m][p],CPLibr1[m],CPLnegNodes[p]); - TSTALLOC(CPLibr2PosPtr[m][p],CPLibr2[m],CPLposNodes[p]); - TSTALLOC(CPLibr1Ibr2Ptr[m][p],CPLibr1[m],CPLibr2[p]); - TSTALLOC(CPLibr2Ibr1Ptr[m][p],CPLibr2[m],CPLibr1[p]); - - } - } - - ReadCpL(here, ckt); - - } - } - - return(OK); + noL = here->dimension; + + here->CPLposNodes = TMALLOC(int, noL); + here->CPLnegNodes = TMALLOC(int, noL); + here->CPLibr1 = TMALLOC(int, noL); + here->CPLibr2 = TMALLOC(int, noL); + + VECTOR_ALLOC(double, here->CPLibr1Ibr1Ptr, noL); + VECTOR_ALLOC(double, here->CPLibr2Ibr2Ptr, noL); + VECTOR_ALLOC(double, here->CPLposIbr1Ptr, noL); + VECTOR_ALLOC(double, here->CPLnegIbr2Ptr, noL); + VECTOR_ALLOC(double, here->CPLposPosPtr, noL); + VECTOR_ALLOC(double, here->CPLnegNegPtr, noL); + VECTOR_ALLOC(double, here->CPLnegPosPtr, noL); + VECTOR_ALLOC(double, here->CPLposNegPtr, noL); + + MATRIX_ALLOC(double, here->CPLibr1PosPtr, noL, noL); + MATRIX_ALLOC(double, here->CPLibr2NegPtr, noL, noL); + MATRIX_ALLOC(double, here->CPLibr1NegPtr, noL, noL); + MATRIX_ALLOC(double, here->CPLibr2PosPtr, noL, noL); + MATRIX_ALLOC(double, here->CPLibr1Ibr2Ptr, noL, noL); + MATRIX_ALLOC(double, here->CPLibr2Ibr1Ptr, noL, noL); + + + branchname = TMALLOC(char*, here->dimension); + if (!here->CPLibr1Given) { + for (m = 0; m < here->dimension; m++) { + branchname[m] = TMALLOC(char, MAX_STRING); + sprintf(branchname[m], "branch1_%d", m); + error = + CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]); + if (error) return (error); + here->CPLibr1[m] = tmp->number; + tfree(branchname[m]); + } + here->CPLibr1Given = 1; + } + tfree(branchname); + branchname = TMALLOC(char*, here->dimension); + if (!here->CPLibr2Given) { + for (m = 0; m < here->dimension; m++) { + branchname[m] = TMALLOC(char, MAX_STRING); + sprintf(branchname[m], "branch2_%d", m); + error = + CKTmkCur(ckt, &tmp, here->CPLname, branchname[m]); + if (error) return (error); + here->CPLibr2[m] = tmp->number; + tfree(branchname[m]); + } + here->CPLibr2Given = 1; + } + tfree(branchname); + + for (m = 0; m < here->dimension; m++) { + for (node = ckt->CKTnodes; node; node = node->next) { + if (strcmp(here->in_node_names[m], + node->name) == 0) { + here->CPLposNodes[m] = node->number; + } + } + } + for (m = 0; m < here->dimension; m++) { + for (node = ckt->CKTnodes; node; node = node->next) { + if (strcmp(here->out_node_names[m], + node->name) == 0) { + here->CPLnegNodes[m] = node->number; + } + } + } + + for (m = 0; m < here->dimension; m++) { + TSTALLOC(CPLibr1Ibr1Ptr[m], CPLibr1[m], CPLibr1[m]); + TSTALLOC(CPLibr2Ibr2Ptr[m], CPLibr2[m], CPLibr2[m]); + TSTALLOC(CPLposIbr1Ptr[m], CPLposNodes[m], CPLibr1[m]); + TSTALLOC(CPLnegIbr2Ptr[m], CPLnegNodes[m], CPLibr2[m]); + TSTALLOC(CPLposPosPtr[m], CPLposNodes[m], CPLposNodes[m]); + TSTALLOC(CPLnegNegPtr[m], CPLnegNodes[m], CPLnegNodes[m]); + TSTALLOC(CPLnegPosPtr[m], CPLnegNodes[m], CPLposNodes[m]); + TSTALLOC(CPLposNegPtr[m], CPLposNodes[m], CPLnegNodes[m]); + + for (p = 0; p < here->dimension; p++) { + + TSTALLOC(CPLibr1PosPtr[m][p], CPLibr1[m], CPLposNodes[p]); + TSTALLOC(CPLibr2NegPtr[m][p], CPLibr2[m], CPLnegNodes[p]); + TSTALLOC(CPLibr1NegPtr[m][p], CPLibr1[m], CPLnegNodes[p]); + TSTALLOC(CPLibr2PosPtr[m][p], CPLibr2[m], CPLposNodes[p]); + TSTALLOC(CPLibr1Ibr2Ptr[m][p], CPLibr1[m], CPLibr2[p]); + TSTALLOC(CPLibr2Ibr1Ptr[m][p], CPLibr2[m], CPLibr1[p]); + + } + } + + ReadCpL(here, ckt); + + } + } + + return(OK); } int -CPLunsetup(GENmodel *inModel, CKTcircuit *ckt) +CPLunsetup(GENmodel* inModel, CKTcircuit* ckt) { - CPLmodel *model; - CPLinstance *here; - int m; - int noL; - - for (model = (CPLmodel *) inModel; model != NULL; - model = CPLnextModel(model)) { - for (here = CPLinstances(model); here != NULL; - here = CPLnextInstance(here)) { - - noL = here->dimension; - - VECTOR_FREE(here->CPLibr1Ibr1Ptr); - VECTOR_FREE(here->CPLibr2Ibr2Ptr); - VECTOR_FREE(here->CPLposIbr1Ptr); - VECTOR_FREE(here->CPLnegIbr2Ptr); - VECTOR_FREE(here->CPLposPosPtr); - VECTOR_FREE(here->CPLnegNegPtr); - VECTOR_FREE(here->CPLnegPosPtr); - VECTOR_FREE(here->CPLposNegPtr); - - - MATRIX_FREE(here->CPLibr1PosPtr, noL, noL); - MATRIX_FREE(here->CPLibr2NegPtr, noL, noL); - MATRIX_FREE(here->CPLibr1NegPtr, noL, noL); - MATRIX_FREE(here->CPLibr2PosPtr, noL, noL); - MATRIX_FREE(here->CPLibr1Ibr2Ptr, noL, noL); - MATRIX_FREE(here->CPLibr2Ibr1Ptr, noL, noL); - - - for (m = 0; m < noL; m++) { - if (here->CPLibr2[m]) { - CKTdltNNum(ckt, here->CPLibr2[m]); - here->CPLibr2[m] = 0; - } - } - - for (m = 0; m < noL; m++) { - if (here->CPLibr1[m]) { - CKTdltNNum(ckt, here->CPLibr1[m]); - here->CPLibr1[m] = 0; - } - } - - tfree(here->CPLposNodes); - tfree(here->CPLnegNodes); - tfree(here->CPLibr1); - tfree(here->CPLibr2); - - /* reset switches */ - here->CPLdcGiven=0; - here->CPLibr1Given = 0; - here->CPLibr2Given = 0; - } - } - mem_delete(); - return OK; + CPLmodel* model; + CPLinstance* here; + int m; + int noL; + + for (model = (CPLmodel*)inModel; model != NULL; + model = CPLnextModel(model)) { + for (here = CPLinstances(model); here != NULL; + here = CPLnextInstance(here)) { + + noL = here->dimension; + + VECTOR_FREE(here->CPLibr1Ibr1Ptr); + VECTOR_FREE(here->CPLibr2Ibr2Ptr); + VECTOR_FREE(here->CPLposIbr1Ptr); + VECTOR_FREE(here->CPLnegIbr2Ptr); + VECTOR_FREE(here->CPLposPosPtr); + VECTOR_FREE(here->CPLnegNegPtr); + VECTOR_FREE(here->CPLnegPosPtr); + VECTOR_FREE(here->CPLposNegPtr); + + + MATRIX_FREE(here->CPLibr1PosPtr, noL, noL); + MATRIX_FREE(here->CPLibr2NegPtr, noL, noL); + MATRIX_FREE(here->CPLibr1NegPtr, noL, noL); + MATRIX_FREE(here->CPLibr2PosPtr, noL, noL); + MATRIX_FREE(here->CPLibr1Ibr2Ptr, noL, noL); + MATRIX_FREE(here->CPLibr2Ibr1Ptr, noL, noL); + + + for (m = 0; m < noL; m++) { + if (here->CPLibr2[m]) { + CKTdltNNum(ckt, here->CPLibr2[m]); + here->CPLibr2[m] = 0; + } + } + + for (m = 0; m < noL; m++) { + if (here->CPLibr1[m]) { + CKTdltNNum(ckt, here->CPLibr1[m]); + here->CPLibr1[m] = 0; + } + } + + tfree(here->CPLposNodes); + tfree(here->CPLnegNodes); + tfree(here->CPLibr1); + tfree(here->CPLibr2); + + /* reset switches */ + here->CPLdcGiven = 0; + here->CPLibr1Given = 0; + here->CPLibr2Given = 0; + } + } + mem_delete(); + return OK; } static int -ReadCpL(CPLinstance *here, CKTcircuit *ckt) +ReadCpL(CPLinstance* here, CKTcircuit* ckt) { - int i, j, noL, counter; - double f; - char *name; - CPLine *c, *c2; - ECPLine *ec; - NODE *nd; - RLINE *lines[MAX_CP_TX_LINES]; - ERLINE *er; - - c = TMALLOC(CPLine, 1); - c2 = TMALLOC(CPLine, 1); - c->vi_head = c->vi_tail = NULL; - noL = c->noL = here->dimension; - here->cplines = c; - here->cplines2 = c2; - - for (i = 0; i < noL; i++) { - ec = TMALLOC(ECPLine, 1); - name = here->in_node_names[i]; - nd = insert_node(name); - ec->link = nd->cplptr; - nd->cplptr = ec; - ec->line = c; - c->in_node[i] = nd; - c2->in_node[i] = nd; - - er = TMALLOC(ERLINE, 1); - er->link = nd->rlptr; - nd->rlptr = er; - er->rl = lines[i] = TMALLOC(RLINE, 1); - er->rl->in_node = nd; - - c->dc1[i] = c->dc2[i] = 0.0; - } - - for (i = 0; i < noL; i++) { - ec = TMALLOC(ECPLine, 1); - name = here->out_node_names[i]; - nd = insert_node(name); - ec->link = nd->cplptr; - nd->cplptr = ec; - ec->line = c; - c->out_node[i] = nd; - c2->out_node[i] = nd; - - er = TMALLOC(ERLINE, 1); - er->link = nd->rlptr; - nd->rlptr = er; - er->rl = lines[i]; - er->rl->out_node = nd; - } - - - counter = 0; - for (i = 0; i < noL; i++) { - for (j = 0; j < noL; j++) { - if (i > j) { - R_m[i][j] = R_m[j][i]; - G_m[i][j] = G_m[j][i]; - C_m[i][j] = C_m[j][i]; - L_m[i][j] = L_m[j][i]; - } else { - f = CPLmodPtr(here)->Rm[counter]; - R_m[i][j] = CPLmodPtr(here)->Rm[counter] = MAX(f, 1.0e-4); - G_m[i][j] = CPLmodPtr(here)->Gm[counter]; - L_m[i][j] = CPLmodPtr(here)->Lm[counter]; - C_m[i][j] = CPLmodPtr(here)->Cm[counter]; - counter++; - } - } - } - if (here->CPLlengthGiven) - length = here->CPLlength; - else length = CPLmodPtr(here)->length; - - for (i = 0; i < noL; i++) - lines[i]->g = 1.0 / (R_m[i][i] * length); - - coupled(noL); - - for (i = 0; i < noL; i++) { - double d, t; - int k; - - c->taul[i] = TAU[i] * 1.0e+12; - for (j = 0; j < noL; j++) { - if (SIV[i][j].C_0 == 0.0) - c->h1t[i][j] = NULL; - else { - c->h1t[i][j] = TMALLOC(TMS, 1); - d = c->h1t[i][j]->aten = SIV[i][j].C_0; - c->h1t[i][j]->ifImg = (int) (SIV[i][j].Poly[6] - 1.0); - /* since originally 2 for img 1 for noimg */ - c->h1t[i][j]->tm[0].c = SIV[i][j].Poly[0] * d; - c->h1t[i][j]->tm[1].c = SIV[i][j].Poly[1] * d; - c->h1t[i][j]->tm[2].c = SIV[i][j].Poly[2] * d; - c->h1t[i][j]->tm[0].x = SIV[i][j].Poly[3]; - c->h1t[i][j]->tm[1].x = SIV[i][j].Poly[4]; - c->h1t[i][j]->tm[2].x = SIV[i][j].Poly[5]; - if (c->h1t[i][j]->ifImg) - c->h1C[i][j] = c->h1t[i][j]->tm[0].c + 2.0 * c->h1t[i][j]->tm[1].c; - else { - t = 0.0; - for (k = 0; k < 3; k++) - t += c->h1t[i][j]->tm[k].c; - c->h1C[i][j] = t; - } - } - - for (k = 0; k < noL; k++) { - if (IWI[i][j].C_0[k] == 0.0) - c->h2t[i][j][k] = NULL; - else { - c->h2t[i][j][k] = TMALLOC(TMS, 1); - d = c->h2t[i][j][k]->aten = IWI[i][j].C_0[k]; - c->h2t[i][j][k]->ifImg = (int) (IWI[i][j].Poly[k][6] - 1.0); - /* since originally 2 for img 1 for noimg */ - c->h2t[i][j][k]->tm[0].c = IWI[i][j].Poly[k][0] * d; - c->h2t[i][j][k]->tm[1].c = IWI[i][j].Poly[k][1] * d; - c->h2t[i][j][k]->tm[2].c = IWI[i][j].Poly[k][2] * d; - c->h2t[i][j][k]->tm[0].x = IWI[i][j].Poly[k][3]; - c->h2t[i][j][k]->tm[1].x = IWI[i][j].Poly[k][4]; - c->h2t[i][j][k]->tm[2].x = IWI[i][j].Poly[k][5]; - if (c->h2t[i][j][k]->ifImg) - c->h2C[i][j][k] = c->h2t[i][j][k]->tm[0].c + 2.0 * - c->h2t[i][j][k]->tm[1].c; - else - c->h2C[i][j][k] = c->h2t[i][j][k]->tm[0].c + - c->h2t[i][j][k]->tm[1].c + - c->h2t[i][j][k]->tm[2].c; - } - if (IWV[i][j].C_0[k] == 0.0) - c->h3t[i][j][k] = NULL; - else { - c->h3t[i][j][k] = TMALLOC(TMS, 1); - d = c->h3t[i][j][k]->aten = IWV[i][j].C_0[k]; - c->h3t[i][j][k]->ifImg = (int) (IWV[i][j].Poly[k][6] - 1.0); - /* since originally 2 for img 1 for noimg */ - c->h3t[i][j][k]->tm[0].c = IWV[i][j].Poly[k][0] * d; - c->h3t[i][j][k]->tm[1].c = IWV[i][j].Poly[k][1] * d; - c->h3t[i][j][k]->tm[2].c = IWV[i][j].Poly[k][2] * d; - c->h3t[i][j][k]->tm[0].x = IWV[i][j].Poly[k][3]; - c->h3t[i][j][k]->tm[1].x = IWV[i][j].Poly[k][4]; - c->h3t[i][j][k]->tm[2].x = IWV[i][j].Poly[k][5]; - if (c->h3t[i][j][k]->ifImg) - c->h3C[i][j][k] = c->h3t[i][j][k]->tm[0].c + 2.0 * - c->h3t[i][j][k]->tm[1].c; - else - c->h3C[i][j][k] = c->h3t[i][j][k]->tm[0].c + - c->h3t[i][j][k]->tm[1].c + - c->h3t[i][j][k]->tm[2].c; - } - } - } - } - - for (i = 0; i < noL; i++) { - if (c->taul[i] < ckt->CKTmaxStep) { - errMsg = TMALLOC(char, strlen(message) + 1); - strcpy(errMsg,message); - return(-1); - } - } - - return(1); + int i, j, noL, counter; + double f; + char* name; + CPLine* c, * c2; + ECPLine* ec; + NODE* nd; + RLINE* lines[MAX_CP_TX_LINES]; + ERLINE* er; + + c = TMALLOC(CPLine, 1); + c2 = TMALLOC(CPLine, 1); + c->vi_head = c->vi_tail = NULL; + noL = c->noL = here->dimension; + here->cplines = c; + here->cplines2 = c2; + + for (i = 0; i < noL; i++) { + ec = TMALLOC(ECPLine, 1); + name = here->in_node_names[i]; + nd = insert_node(name); + ec->link = nd->cplptr; + nd->cplptr = ec; + ec->line = c; + c->in_node[i] = nd; + c2->in_node[i] = nd; + + er = TMALLOC(ERLINE, 1); + er->link = nd->rlptr; + nd->rlptr = er; + er->rl = lines[i] = TMALLOC(RLINE, 1); + er->rl->in_node = nd; + + c->dc1[i] = c->dc2[i] = 0.0; + } + + for (i = 0; i < noL; i++) { + ec = TMALLOC(ECPLine, 1); + name = here->out_node_names[i]; + nd = insert_node(name); + ec->link = nd->cplptr; + nd->cplptr = ec; + ec->line = c; + c->out_node[i] = nd; + c2->out_node[i] = nd; + + er = TMALLOC(ERLINE, 1); + er->link = nd->rlptr; + nd->rlptr = er; + er->rl = lines[i]; + er->rl->out_node = nd; + } + + + counter = 0; + for (i = 0; i < noL; i++) { + for (j = 0; j < noL; j++) { + if (i > j) { + R_m[i][j] = R_m[j][i]; + G_m[i][j] = G_m[j][i]; + C_m[i][j] = C_m[j][i]; + L_m[i][j] = L_m[j][i]; + } + else { + f = CPLmodPtr(here)->Rm[counter]; + R_m[i][j] = CPLmodPtr(here)->Rm[counter] = MAX(f, 1.0e-4); + G_m[i][j] = CPLmodPtr(here)->Gm[counter]; + L_m[i][j] = CPLmodPtr(here)->Lm[counter]; + C_m[i][j] = CPLmodPtr(here)->Cm[counter]; + counter++; + } + } + } + if (here->CPLlengthGiven) + length = here->CPLlength; + else length = CPLmodPtr(here)->length; + + for (i = 0; i < noL; i++) + lines[i]->g = 1.0 / (R_m[i][i] * length); + + coupled(noL); + + for (i = 0; i < noL; i++) { + double d, t; + int k; + + c->taul[i] = TAU[i] * 1.0e+12; + for (j = 0; j < noL; j++) { + if (SIV[i][j].C_0 == 0.0) + c->h1t[i][j] = NULL; + else { + c->h1t[i][j] = TMALLOC(TMS, 1); + d = c->h1t[i][j]->aten = SIV[i][j].C_0; + c->h1t[i][j]->ifImg = (int)(SIV[i][j].Poly[6] - 1.0); + /* since originally 2 for img 1 for noimg */ + c->h1t[i][j]->tm[0].c = SIV[i][j].Poly[0] * d; + c->h1t[i][j]->tm[1].c = SIV[i][j].Poly[1] * d; + c->h1t[i][j]->tm[2].c = SIV[i][j].Poly[2] * d; + c->h1t[i][j]->tm[0].x = SIV[i][j].Poly[3]; + c->h1t[i][j]->tm[1].x = SIV[i][j].Poly[4]; + c->h1t[i][j]->tm[2].x = SIV[i][j].Poly[5]; + if (c->h1t[i][j]->ifImg) + c->h1C[i][j] = c->h1t[i][j]->tm[0].c + 2.0 * c->h1t[i][j]->tm[1].c; + else { + t = 0.0; + for (k = 0; k < 3; k++) + t += c->h1t[i][j]->tm[k].c; + c->h1C[i][j] = t; + } + } + + for (k = 0; k < noL; k++) { + if (IWI[i][j].C_0[k] == 0.0) + c->h2t[i][j][k] = NULL; + else { + c->h2t[i][j][k] = TMALLOC(TMS, 1); + d = c->h2t[i][j][k]->aten = IWI[i][j].C_0[k]; + c->h2t[i][j][k]->ifImg = (int)(IWI[i][j].Poly[k][6] - 1.0); + /* since originally 2 for img 1 for noimg */ + c->h2t[i][j][k]->tm[0].c = IWI[i][j].Poly[k][0] * d; + c->h2t[i][j][k]->tm[1].c = IWI[i][j].Poly[k][1] * d; + c->h2t[i][j][k]->tm[2].c = IWI[i][j].Poly[k][2] * d; + c->h2t[i][j][k]->tm[0].x = IWI[i][j].Poly[k][3]; + c->h2t[i][j][k]->tm[1].x = IWI[i][j].Poly[k][4]; + c->h2t[i][j][k]->tm[2].x = IWI[i][j].Poly[k][5]; + if (c->h2t[i][j][k]->ifImg) + c->h2C[i][j][k] = c->h2t[i][j][k]->tm[0].c + 2.0 * + c->h2t[i][j][k]->tm[1].c; + else + c->h2C[i][j][k] = c->h2t[i][j][k]->tm[0].c + + c->h2t[i][j][k]->tm[1].c + + c->h2t[i][j][k]->tm[2].c; + } + if (IWV[i][j].C_0[k] == 0.0) + c->h3t[i][j][k] = NULL; + else { + c->h3t[i][j][k] = TMALLOC(TMS, 1); + d = c->h3t[i][j][k]->aten = IWV[i][j].C_0[k]; + c->h3t[i][j][k]->ifImg = (int)(IWV[i][j].Poly[k][6] - 1.0); + /* since originally 2 for img 1 for noimg */ + c->h3t[i][j][k]->tm[0].c = IWV[i][j].Poly[k][0] * d; + c->h3t[i][j][k]->tm[1].c = IWV[i][j].Poly[k][1] * d; + c->h3t[i][j][k]->tm[2].c = IWV[i][j].Poly[k][2] * d; + c->h3t[i][j][k]->tm[0].x = IWV[i][j].Poly[k][3]; + c->h3t[i][j][k]->tm[1].x = IWV[i][j].Poly[k][4]; + c->h3t[i][j][k]->tm[2].x = IWV[i][j].Poly[k][5]; + if (c->h3t[i][j][k]->ifImg) + c->h3C[i][j][k] = c->h3t[i][j][k]->tm[0].c + 2.0 * + c->h3t[i][j][k]->tm[1].c; + else + c->h3C[i][j][k] = c->h3t[i][j][k]->tm[0].c + + c->h3t[i][j][k]->tm[1].c + + c->h3t[i][j][k]->tm[2].c; + } + } + } + } + + for (i = 0; i < noL; i++) { + if (c->taul[i] < ckt->CKTmaxStep) { + errMsg = TMALLOC(char, strlen(message) + 1); + strcpy(errMsg, message); + return(-1); + } + } + + return(1); } /**************************************************************** - misc.c Miscellaneous procedures for simulation of - coupled transmission lines. + misc.c Miscellaneous procedures for simulation of + coupled transmission lines. ****************************************************************/ static void new_memory(int dim, int deg, int deg_o) { - int i, j; + int i, j; - NG_IGNORE(deg); + NG_IGNORE(deg); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - SiSv_1[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + SiSv_1[i][j] = (double*)calloc((size_t)(deg_o + 1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sip[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Sip[i][j] = (double*)calloc((size_t)(deg_o + 1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si_1p[i][j] = (double*)calloc((size_t)(deg_o + 1), sizeof(double)); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sv_1p[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Sv_1p[i][j] = (double*)calloc((size_t)(deg_o + 1), sizeof(double)); - for (i = 0; i < dim; i++) - W[i] = (double *) calloc(MAX_DEG, sizeof(double)); + for (i = 0; i < dim; i++) + W[i] = (double*)calloc(MAX_DEG, sizeof(double)); } /*** ***/ -/**************************************************************** - match Create a polynomial matching given data points - ****************************************************************/ + /**************************************************************** + match Create a polynomial matching given data points + ****************************************************************/ -static double * +static double* vector(int nl, int nh) { - double *v = TMALLOC(double, nh - nl + 1); + double* v = TMALLOC(double, nh - nl + 1); - if (!v) { - fprintf(stderr, "Memory Allocation Error by tmalloc in vector().\n"); - fprintf(stderr, "...now exiting to system ...\n"); - controlled_exit(EXIT_FAILURE); - } + if (!v) { + fprintf(stderr, "Memory Allocation Error by tmalloc in vector().\n"); + fprintf(stderr, "...now exiting to system ...\n"); + controlled_exit(EXIT_FAILURE); + } - return v - nl; + return v - nl; } static void -free_vector(double *v, int nl, int nh) +free_vector(double* v, int nl, int nh) { - NG_IGNORE(nh); - double *freev = v + nl; - tfree(freev); + NG_IGNORE(nh); + double* freev = v + nl; + tfree(freev); } static void -polint(double *xa, double *ya, int n, double x, double *y, double *dy) +polint(double* xa, double* ya, int n, double x, double* y, double* dy) /* Given arrays xa[1..n] and ya[1..n], and given a value x, this routine returns a value y, and an error estimate dy. If P(x) is the @@ -601,199 +602,199 @@ polint(double *xa, double *ya, int n, double x, double *y, double *dy) value y = P(x) */ { - int i, m, ns = 1; - double den, dif, dift, ho, hp, w; - double *c, *d; - - dif = ABS(x - xa[1]); - c = vector(1, n); - d = vector(1, n); - for (i = 1; i <= n; i++) { - if ((dift = ABS(x - xa[i])) < dif) { - ns = i; - dif = dift; - } - c[i] = ya[i]; - d[i] = ya[i]; - } - *y = ya[ns--]; - for (m = 1; m < n; m++) { - for (i = 1; i <= n-m; i++) { - ho = xa[i]-x; - hp = xa[i+m]-x; - w = c[i+1]-d[i]; - if ((den=ho-hp) == 0.0) { - fprintf(stderr, "(Error) in routine POLINT\n"); - fprintf(stderr, "...now exiting to system ...\n"); - controlled_exit(EXIT_FAILURE); - } - den = w/den; - d[i] = hp * den; - c[i] = ho * den; - } - *y += (*dy = (2*ns < (n-m) ? c[ns+1] : d[ns--])); - } - free_vector(d, 1, n); - free_vector(c, 1, n); + int i, m, ns = 1; + double den, dif, dift, ho, hp, w; + double* c, * d; + + dif = ABS(x - xa[1]); + c = vector(1, n); + d = vector(1, n); + for (i = 1; i <= n; i++) { + if ((dift = ABS(x - xa[i])) < dif) { + ns = i; + dif = dift; + } + c[i] = ya[i]; + d[i] = ya[i]; + } + *y = ya[ns--]; + for (m = 1; m < n; m++) { + for (i = 1; i <= n - m; i++) { + ho = xa[i] - x; + hp = xa[i + m] - x; + w = c[i + 1] - d[i]; + if ((den = ho - hp) == 0.0) { + fprintf(stderr, "(Error) in routine POLINT\n"); + fprintf(stderr, "...now exiting to system ...\n"); + controlled_exit(EXIT_FAILURE); + } + den = w / den; + d[i] = hp * den; + c[i] = ho * den; + } + *y += (*dy = (2 * ns < (n - m) ? c[ns + 1] : d[ns--])); + } + free_vector(d, 1, n); + free_vector(c, 1, n); } static int -match(int n, double *cof, double *xa, double *ya) +match(int n, double* cof, double* xa, double* ya) /* Given arrays xa[0..n] and ya[0..n] containing a tabulated function ya = f(xa), this routine returns an array of coefficients cof[0..n], such that ya[i] = sum_j {cof[j]*xa[i]**j}. */ { - int k, j, i; - double xmin, dy, *x, *y, *xx; - - x = vector(0, n); - y = vector(0, n); - xx = vector(0, n); - for (j = 0; j <= n; j++) { - x[j] = xa[j]; - xx[j] = y[j] = ya[j]; - } - for (j = 0; j <= n; j++) { - polint(x-1, y-1, n+1-j, 0.0, &cof[j], &dy); - xmin = 1.0e38; - k = -1; - for (i = 0; i <= n-j; i++) { - if (ABS(x[i]) < xmin) { - xmin = ABS(x[i]); - k = i; - } - if (x[i]) y[i] = (y[i] - cof[j]) / x[i]; - } - for (i = k+1; i <= n-j; i++) { - y[i-1] = y[i]; - x[i-1] = x[i]; - } - } - free_vector(y, 0, n); - free_vector(x, 0, n); - free_vector(xx, 0, n); - - /**** check ****/ - /* - for (i = 0; i <= n; i++) { - xmin = xa[i]; - dy = cof[0]; - for (j = 1; j <= n; j++) { - dy += xmin * cof[j]; - xmin *= xa[i]; - } - printf("*** real x = %e y = %e\n", xa[i], xx[i]); - printf("*** calculated y = %e\n", dy); - printf("*** error = %e \% \n", (dy-xx[i])/xx[i]); - } - */ - return 0; + int k, j, i; + double xmin, dy, * x, * y, * xx; + + x = vector(0, n); + y = vector(0, n); + xx = vector(0, n); + for (j = 0; j <= n; j++) { + x[j] = xa[j]; + xx[j] = y[j] = ya[j]; + } + for (j = 0; j <= n; j++) { + polint(x - 1, y - 1, n + 1 - j, 0.0, &cof[j], &dy); + xmin = 1.0e38; + k = -1; + for (i = 0; i <= n - j; i++) { + if (ABS(x[i]) < xmin) { + xmin = ABS(x[i]); + k = i; + } + if (x[i]) y[i] = (y[i] - cof[j]) / x[i]; + } + for (i = k + 1; i <= n - j; i++) { + y[i - 1] = y[i]; + x[i - 1] = x[i]; + } + } + free_vector(y, 0, n); + free_vector(x, 0, n); + free_vector(xx, 0, n); + + /**** check ****/ + /* + for (i = 0; i <= n; i++) { + xmin = xa[i]; + dy = cof[0]; + for (j = 1; j <= n; j++) { + dy += xmin * cof[j]; + xmin *= xa[i]; + } + printf("*** real x = %e y = %e\n", xa[i], xx[i]); + printf("*** calculated y = %e\n", dy); + printf("*** error = %e \% \n", (dy-xx[i])/xx[i]); + } + */ + return 0; } /*** ***/ -/*** -static int -match_x(int dim, double *Alfa, double *X, double *Y) -{ - int i, j; - double f; - double scale; - - **** check **** - double xx[16]; - for (i = 0; i <= dim; i++) - xx[i] = Y[i]; - - if (Y[1] == Y[0]) - scale = 1.0; - else - scale = X[1] / (Y[1] - Y[0]); - for (i = 0; i < dim; i++) { - f = X[i+1]; - for (j = dim-1; j >= 0; j--) { - A[i][j] = f; - f *= X[i+1]; - } - A[i][dim] = (Y[i+1] - Y[0])*scale; - } - Gaussian_Elimination2(dim, 1); - Alfa[0] = Y[0]; - for (i = 1; i <= dim; i++) - Alfa[i] = A[dim-i][dim] / scale; - - **** check **** - * - for (i = 0; i <= dim; i++) { - f = X[i]; - scale = Alfa[0]; - for (j = 1; j <= dim; j++) { - scale += f * Alfa[j]; - f *= X[i]; - } - printf("*** real x = %e y = %e\n", X[i], xx[i]); - printf("*** calculated y = %e\n", scale); - printf("*** error = %e \% \n", (scale-xx[i])/xx[i]); - } - * - - return(1); -} -***/ -/*** + /*** + static int + match_x(int dim, double *Alfa, double *X, double *Y) + { + int i, j; + double f; + double scale; + + **** check **** + double xx[16]; + for (i = 0; i <= dim; i++) + xx[i] = Y[i]; + + if (Y[1] == Y[0]) + scale = 1.0; + else + scale = X[1] / (Y[1] - Y[0]); + for (i = 0; i < dim; i++) { + f = X[i+1]; + for (j = dim-1; j >= 0; j--) { + A[i][j] = f; + f *= X[i+1]; + } + A[i][dim] = (Y[i+1] - Y[0])*scale; + } + Gaussian_Elimination2(dim, 1); + Alfa[0] = Y[0]; + for (i = 1; i <= dim; i++) + Alfa[i] = A[dim-i][dim] / scale; + + **** check **** + * + for (i = 0; i <= dim; i++) { + f = X[i]; + scale = Alfa[0]; + for (j = 1; j <= dim; j++) { + scale += f * Alfa[j]; + f *= X[i]; + } + printf("*** real x = %e y = %e\n", X[i], xx[i]); + printf("*** calculated y = %e\n", scale); + printf("*** error = %e \% \n", (scale-xx[i])/xx[i]); + } + * + + return(1); + } ***/ + /*** + ***/ static int Gaussian_Elimination2(int dims, int type) /* type = 1 : to solve a linear system - -1 : to inverse a matrix */ + -1 : to inverse a matrix */ { - int i, j, k, dim; - double f; - double max; - int imax; - - if (type == -1) - dim = 2 * dims; - else - dim = dims; - - for (i = 0; i < dims; i++) { - imax = i; - max = ABS(A[i][i]); - for (j = i+1; j < dim; j++) - if (ABS(A[j][i]) > max) { - imax = j; - max = ABS(A[j][i]); - } - if (max < epsilon) { - fprintf(stderr, " can not choose a pivot (misc)\n"); - controlled_exit(EXIT_FAILURE); - } - if (imax != i) - for (k = i; k <= dim; k++) { - SWAP(double, A[i][k], A[imax][k]); - } - - f = 1.0 / A[i][i]; - A[i][i] = 1.0; - - for (j = i+1; j <= dim; j++) - A[i][j] *= f; - - for (j = 0; j < dims ; j++) { - if (i == j) - continue; - f = A[j][i]; - A[j][i] = 0.0; - for (k = i+1; k <= dim; k++) - A[j][k] -= f * A[i][k]; - } - } - - return(1); + int i, j, k, dim; + double f; + double max; + int imax; + + if (type == -1) + dim = 2 * dims; + else + dim = dims; + + for (i = 0; i < dims; i++) { + imax = i; + max = ABS(A[i][i]); + for (j = i + 1; j < dim; j++) + if (ABS(A[j][i]) > max) { + imax = j; + max = ABS(A[j][i]); + } + if (max < epsilon) { + fprintf(stderr, " can not choose a pivot (misc)\n"); + controlled_exit(EXIT_FAILURE); + } + if (imax != i) + for (k = i; k <= dim; k++) { + SWAP(double, A[i][k], A[imax][k]); + } + + f = 1.0 / A[i][i]; + A[i][i] = 1.0; + + for (j = i + 1; j <= dim; j++) + A[i][j] *= f; + + for (j = 0; j < dims; j++) { + if (i == j) + continue; + f = A[j][i]; + A[j][i] = 0.0; + for (k = i + 1; k <= dim; k++) + A[j][k] -= f * A[i][k]; + } + } + + return(1); } /*** @@ -804,35 +805,35 @@ eval_Si_Si_1(int dim, double y) int i, j, k; for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Si_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - if (k == j) - Si_1[i][j] += Sv_1[i][k] * - (y * R_m[k][j] + Scaling_F * L_m[k][j]); - else - Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F; - / - Si_1[i][j] *= Scaling_F; - / - } + for (j = 0; j < dim; j++) { + Si_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + if (k == j) + Si_1[i][j] += Sv_1[i][k] * + (y * R_m[k][j] + Scaling_F * L_m[k][j]); + else + Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F; + / + Si_1[i][j] *= Scaling_F; + / + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si_1[i][j] /= sqrt((double) D[i]); + for (j = 0; j < dim; j++) + Si_1[i][j] /= sqrt((double) D[i]); for (i = 0; i < dim; i++) { - for (j = 0; j < dim; j++) - A[i][j] = Si_1[i][j]; - for (j = dim; j < 2* dim; j++) - A[i][j] = 0.0; - A[i][i+dim] = 1.0; + for (j = 0; j < dim; j++) + A[i][j] = Si_1[i][j]; + for (j = dim; j < 2* dim; j++) + A[i][j] = 0.0; + A[i][i+dim] = 1.0; } Gaussian_Elimination2(dim, -1); for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si[i][j] = A[i][j+dim]; + for (j = 0; j < dim; j++) + Si[i][j] = A[i][j+dim]; } ***/ @@ -840,36 +841,36 @@ eval_Si_Si_1(int dim, double y) static void eval_Si_Si_1(int dim, double y) { - int i, j, k; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Si_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Si_1[i][j] += Sv_1[i][k] * (y * R_m[k][j] + Scaling_F * L_m[k][j]); - /* - else - Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F; - Si_1[i][j] *= Scaling_F; - */ - } - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si_1[i][j] /= sqrt(D[i]); - - for (i = 0; i < dim; i++) { - for (j = 0; j < dim; j++) - A[i][j] = Si_1[i][j]; - for (j = dim; j < 2* dim; j++) - A[i][j] = 0.0; - A[i][i+dim] = 1.0; - } - Gaussian_Elimination2(dim, -1); - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Si[i][j] = A[i][j+dim]; + int i, j, k; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Si_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Si_1[i][j] += Sv_1[i][k] * (y * R_m[k][j] + Scaling_F * L_m[k][j]); + /* + else + Si_1[i][j] += Sv_1[i][k] * L_m[k][j] * Scaling_F; + Si_1[i][j] *= Scaling_F; + */ + } + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si_1[i][j] /= sqrt(D[i]); + + for (i = 0; i < dim; i++) { + for (j = 0; j < dim; j++) + A[i][j] = Si_1[i][j]; + for (j = dim; j < 2 * dim; j++) + A[i][j] = 0.0; + A[i][i + dim] = 1.0; + } + Gaussian_Elimination2(dim, -1); + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Si[i][j] = A[i][j + dim]; } /*** @@ -881,88 +882,88 @@ loop_ZY(int dim, double y) double fmin, fmin1; for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - if (i == j) - ZY[i][j] = Scaling_F * C_m[i][i] + G_m[i] * y; - else - ZY[i][j] = Scaling_F * C_m[i][j]; + for (j = 0; j < dim; j++) + if (i == j) + ZY[i][j] = Scaling_F * C_m[i][i] + G_m[i] * y; + else + ZY[i][j] = Scaling_F * C_m[i][j]; diag(dim); fmin = D[0]; for (i = 1; i < dim; i++) - if (D[i] < fmin) - fmin = D[i]; + if (D[i] < fmin) + fmin = D[i]; if (fmin < 0) { - fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n"); - exit(0); + fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n"); + exit(0); } else { - fmin = sqrt(fmin); - fmin1 = 1 / fmin; + fmin = sqrt(fmin); + fmin1 = 1 / fmin; } for (i = 0; i < dim; i++) - D[i] = sqrt((double) D[i]); + D[i] = sqrt((double) D[i]); for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Y5[i][j] = D[i] * Sv[j][i]; - Y5_1[i][j] = Sv[j][i] / D[i]; - } + for (j = 0; j < dim; j++) { + Y5[i][j] = D[i] * Sv[j][i]; + Y5_1[i][j] = Sv[j][i] / D[i]; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[i][k] * Y5[k][j]; - } + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[i][k] * Y5[k][j]; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Y5[i][j] = Sv_1[i][j]; + for (j = 0; j < dim; j++) + Y5[i][j] = Sv_1[i][j]; for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[i][k] * Y5_1[k][j]; - } + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[i][k] * Y5_1[k][j]; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Y5_1[i][j] = Sv_1[i][j]; + for (j = 0; j < dim; j++) + Y5_1[i][j] = Sv_1[i][j]; for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - ZY[i][j] = 0.0; - for (k = 0; k < dim; k++) - if (k == i) - ZY[i][j] += (Scaling_F * L_m[i][i] + R_m[i] * y) * - Y5[k][j]; - else - ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F; - } + for (j = 0; j < dim; j++) { + ZY[i][j] = 0.0; + for (k = 0; k < dim; k++) + if (k == i) + ZY[i][j] += (Scaling_F * L_m[i][i] + R_m[i] * y) * + Y5[k][j]; + else + ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Y5[i][k] * ZY[k][j]; - } + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Y5[i][k] * ZY[k][j]; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - ZY[i][j] = Sv_1[i][j]; + for (j = 0; j < dim; j++) + ZY[i][j] = Sv_1[i][j]; diag(dim); for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[k][i] * Y5[k][j]; - Sv_1[i][j] *= fmin1; - } + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[k][i] * Y5[k][j]; + Sv_1[i][j] *= fmin1; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - ZY[i][j] = 0.0; - for (k = 0; k < dim; k++) - ZY[i][j] += Y5_1[i][k] * Sv[k][j]; - ZY[i][j] *= fmin; - } + for (j = 0; j < dim; j++) { + ZY[i][j] = 0.0; + for (k = 0; k < dim; k++) + ZY[i][j] += Y5_1[i][k] * Sv[k][j]; + ZY[i][j] *= fmin; + } for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sv[i][j] = ZY[i][j]; + for (j = 0; j < dim; j++) + Sv[i][j] = ZY[i][j]; } ***/ @@ -970,93 +971,94 @@ loop_ZY(int dim, double y) static void loop_ZY(int dim, double y) { - int i, j, k; - double fmin, fmin1=0.0; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - ZY[i][j] = Scaling_F * C_m[i][j] + G_m[i][j] * y; - /* - else - ZY[i][j] = Scaling_F * C_m[i][j]; - */ - diag(dim); - fmin = D[0]; - for (i = 1; i < dim; i++) - if (D[i] < fmin) - fmin = D[i]; - if (fmin < 0) { - fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n"); - controlled_exit(EXIT_FAILURE); - } else { - fmin = sqrt(fmin); - fmin1 = 1 / fmin; - } - for (i = 0; i < dim; i++) - D[i] = sqrt(D[i]); - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Y5[i][j] = D[i] * Sv[j][i]; - Y5_1[i][j] = Sv[j][i] / D[i]; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[i][k] * Y5[k][j]; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Y5[i][j] = Sv_1[i][j]; - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[i][k] * Y5_1[k][j]; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Y5_1[i][j] = Sv_1[i][j]; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - ZY[i][j] = 0.0; - for (k = 0; k < dim; k++) - ZY[i][j] += (Scaling_F * L_m[i][k] + R_m[i][k] * y) * Y5[k][j]; - /* - else - ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F; - */ - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Y5[i][k] * ZY[k][j]; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - ZY[i][j] = Sv_1[i][j]; - - diag(dim); - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - Sv_1[i][j] = 0.0; - for (k = 0; k < dim; k++) - Sv_1[i][j] += Sv[k][i] * Y5[k][j]; - Sv_1[i][j] *= fmin1; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - ZY[i][j] = 0.0; - for (k = 0; k < dim; k++) - ZY[i][j] += Y5_1[i][k] * Sv[k][j]; - ZY[i][j] *= fmin; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - Sv[i][j] = ZY[i][j]; + int i, j, k; + double fmin, fmin1 = 0.0; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + ZY[i][j] = Scaling_F * C_m[i][j] + G_m[i][j] * y; + /* + else + ZY[i][j] = Scaling_F * C_m[i][j]; + */ + diag(dim); + fmin = D[0]; + for (i = 1; i < dim; i++) + if (D[i] < fmin) + fmin = D[i]; + if (fmin < 0) { + fprintf(stderr, "(Error) The capacitance matrix of the multiconductor system is not positive definite.\n"); + controlled_exit(EXIT_FAILURE); + } + else { + fmin = sqrt(fmin); + fmin1 = 1 / fmin; + } + for (i = 0; i < dim; i++) + D[i] = sqrt(D[i]); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Y5[i][j] = D[i] * Sv[j][i]; + Y5_1[i][j] = Sv[j][i] / D[i]; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[i][k] * Y5[k][j]; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Y5[i][j] = Sv_1[i][j]; + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[i][k] * Y5_1[k][j]; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Y5_1[i][j] = Sv_1[i][j]; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + ZY[i][j] = 0.0; + for (k = 0; k < dim; k++) + ZY[i][j] += (Scaling_F * L_m[i][k] + R_m[i][k] * y) * Y5[k][j]; + /* + else + ZY[i][j] += L_m[i][k] * Y5[k][j] * Scaling_F; + */ + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Y5[i][k] * ZY[k][j]; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + ZY[i][j] = Sv_1[i][j]; + + diag(dim); + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + Sv_1[i][j] = 0.0; + for (k = 0; k < dim; k++) + Sv_1[i][j] += Sv[k][i] * Y5[k][j]; + Sv_1[i][j] *= fmin1; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + ZY[i][j] = 0.0; + for (k = 0; k < dim; k++) + ZY[i][j] += Y5_1[i][k] * Sv[k][j]; + ZY[i][j] *= fmin; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + Sv[i][j] = ZY[i][j]; } @@ -1066,56 +1068,56 @@ loop_ZY(int dim, double y) static void poly_matrix( - double *A_in[MAX_DIM][MAX_DIM], - int dim, int deg) + double* A_in[MAX_DIM][MAX_DIM], + int dim, int deg) { - int i, j; + int i, j; - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - match(deg, A_in[i][j], frequency, A_in[i][j]); + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + match(deg, A_in[i][j], frequency, A_in[i][j]); } /*** ***/ -/*** -static int -checkW(double *W, double d) -{ - double f, y; - float y1; - int k; - - printf("(W)y ="); - scanf("%f", &y1); - - f = W[0]; - y = y1; - f += y * W[1]; - for (k = 2; k < 6; k++) { - y *= y1; - f += y * W[k]; - } - printf("W[i]= %e\n ", f*exp((double)-d/y1)); - - return(1); -} -***/ -/*** + /*** + static int + checkW(double *W, double d) + { + double f, y; + float y1; + int k; + + printf("(W)y ="); + scanf("%f", &y1); + + f = W[0]; + y = y1; + f += y * W[1]; + for (k = 2; k < 6; k++) { + y *= y1; + f += y * W[k]; + } + printf("W[i]= %e\n ", f*exp((double)-d/y1)); + + return(1); + } ***/ + /*** + ***/ static void poly_W(int dim, int deg) { - int i; - - for (i = 0; i < dim; i++) { - match(deg, W[i], frequency, W[i]); - TAU[i] = approx_mode(W[i], W[i], length); - /* - checkW(W[i], TAU[i]); - */ - } + int i; + + for (i = 0; i < dim; i++) { + match(deg, W[i], frequency, W[i]); + TAU[i] = approx_mode(W[i], W[i], length); + /* + checkW(W[i], TAU[i]); + */ + } } /*** @@ -1124,36 +1126,36 @@ poly_W(int dim, int deg) static void eval_frequency(int dim, int deg_o) { - int i; - double min; + int i; + double min; - min = D[0]; + min = D[0]; - for (i = 1; i < dim; i++) - if (D[i] < min) { - min = D[i]; - } + for (i = 1; i < dim; i++) + if (D[i] < min) { + min = D[i]; + } - if (min <= 0) { - fprintf(stderr, "A mode frequency is not positive. Abort!\n"); - controlled_exit(EXIT_FAILURE); - } + if (min <= 0) { + fprintf(stderr, "A mode frequency is not positive. Abort!\n"); + controlled_exit(EXIT_FAILURE); + } - Scaling_F2 = 1.0 / min; - Scaling_F = sqrt(Scaling_F2); - min = length * 8.0; - /* - min *= 1.0e18; - min = sqrt(min)*1.0e-9*length/8.0; - */ + Scaling_F2 = 1.0 / min; + Scaling_F = sqrt(Scaling_F2); + min = length * 8.0; + /* + min *= 1.0e18; + min = sqrt(min)*1.0e-9*length/8.0; + */ - frequency[0] = 0.0; + frequency[0] = 0.0; - for (i = 1; i <= deg_o; i++) - frequency[i] = frequency[i-1] + min; + for (i = 1; i <= deg_o; i++) + frequency[i] = frequency[i - 1] + min; - for (i = 0; i < dim; i++) - D[i] *= Scaling_F2; + for (i = 0; i < dim; i++) + D[i] *= Scaling_F2; } /*** @@ -1162,20 +1164,20 @@ eval_frequency(int dim, int deg_o) static void store(int dim, int ind) { - int i, j; - - for (i = 0; i < dim; i++) { - for (j = 0; j < dim; j++) { - /* store_Sip */ - Sip[i][j][ind] = Si[i][j]; - /* store_Si_1p */ - Si_1p[i][j][ind] = Si_1[i][j]; - /* store_Sv_1p */ - Sv_1p[i][j][ind] = Sv_1[i][j]; - } - /* store_W */ - W[i][ind] = D[i]; - } + int i, j; + + for (i = 0; i < dim; i++) { + for (j = 0; j < dim; j++) { + /* store_Sip */ + Sip[i][j][ind] = Si[i][j]; + /* store_Si_1p */ + Si_1p[i][j][ind] = Si_1[i][j]; + /* store_Sv_1p */ + Sv_1p[i][j][ind] = Sv_1[i][j]; + } + /* store_W */ + W[i][ind] = D[i]; + } } /*** @@ -1184,144 +1186,144 @@ store(int dim, int ind) static void store_SiSv_1(int dim, int ind) { - int i, j, k; - double temp; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - temp = 0.0; - for (k = 0; k < dim; k++) - temp += Si[i][k] * Sv_1[k][j]; - SiSv_1[i][j][ind] = temp; - } + int i, j, k; + double temp; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + temp = 0.0; + for (k = 0; k < dim; k++) + temp += Si[i][k] * Sv_1[k][j]; + SiSv_1[i][j][ind] = temp; + } } /*** ***/ -/*** -static int -check(Sip, Si_1p, Sv_1p, SiSv_1p) - double *Sip[MAX_DIM][MAX_DIM]; - double *Si_1p[MAX_DIM][MAX_DIM]; - double *Sv_1p[MAX_DIM][MAX_DIM]; - double *SiSv_1p[MAX_DIM][MAX_DIM]; -{ - double f, y; - float y1; - int i, j, k; - - printf("y ="); - scanf("%f", &y1); - - printf("\n"); - printf("Si =\n"); - for (i = 0; i < 4; i++) { - for (j = 0; j < 4; j++) { - f = Sip[i][j][0]; - y = y1; - f += y * Sip[i][j][1]; - for (k = 2; k < 8; k++) { - y *= y1; - f += y * Sip[i][j][k]; - } - printf("%e ", f); - } - printf("\n"); - } - printf("\n"); - printf("Si_1 =\n"); - for (i = 0; i < 4; i++) { - for (j = 0; j < 4; j++) { - f = Si_1p[i][j][0]; - y = y1; - f += y * Si_1p[i][j][1]; - for (k = 2; k < 8; k++) { - y *= y1; - f += y * Si_1p[i][j][k]; - } - printf("%e ", f); - } - printf("\n"); - } - printf("\n"); - printf("Sv_1 =\n"); - for (i = 0; i < 4; i++) { - for (j = 0; j < 4; j++) { - f = Sv_1p[i][j][0]; - y = y1; - f += y * Sv_1p[i][j][1]; - for (k = 2; k < 8; k++) { - y *= y1; - f += y * Sv_1p[i][j][k]; - } - printf("%e ", f); - } - printf("\n"); - } - printf("\n"); - printf("SiSv_1 =\n"); - for (i = 0; i < 4; i++) { - for (j = 0; j < 4; j++) { - f = SiSv_1p[i][j][0]; - y = y1; - f += y * SiSv_1p[i][j][1]; - for (k = 2; k < 8; k++) { - y *= y1; - f += y * SiSv_1p[i][j][k]; - } - printf("%e ", f); - } - printf("\n"); - } - return(1); -} -***/ -/*** + /*** + static int + check(Sip, Si_1p, Sv_1p, SiSv_1p) + double *Sip[MAX_DIM][MAX_DIM]; + double *Si_1p[MAX_DIM][MAX_DIM]; + double *Sv_1p[MAX_DIM][MAX_DIM]; + double *SiSv_1p[MAX_DIM][MAX_DIM]; + { + double f, y; + float y1; + int i, j, k; + + printf("y ="); + scanf("%f", &y1); + + printf("\n"); + printf("Si =\n"); + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + f = Sip[i][j][0]; + y = y1; + f += y * Sip[i][j][1]; + for (k = 2; k < 8; k++) { + y *= y1; + f += y * Sip[i][j][k]; + } + printf("%e ", f); + } + printf("\n"); + } + printf("\n"); + printf("Si_1 =\n"); + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + f = Si_1p[i][j][0]; + y = y1; + f += y * Si_1p[i][j][1]; + for (k = 2; k < 8; k++) { + y *= y1; + f += y * Si_1p[i][j][k]; + } + printf("%e ", f); + } + printf("\n"); + } + printf("\n"); + printf("Sv_1 =\n"); + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + f = Sv_1p[i][j][0]; + y = y1; + f += y * Sv_1p[i][j][1]; + for (k = 2; k < 8; k++) { + y *= y1; + f += y * Sv_1p[i][j][k]; + } + printf("%e ", f); + } + printf("\n"); + } + printf("\n"); + printf("SiSv_1 =\n"); + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + f = SiSv_1p[i][j][0]; + y = y1; + f += y * SiSv_1p[i][j][1]; + for (k = 2; k < 8; k++) { + y *= y1; + f += y * SiSv_1p[i][j][k]; + } + printf("%e ", f); + } + printf("\n"); + } + return(1); + } ***/ + /*** + ***/ static int coupled(int dim) { - int deg, deg_o; - int i; - - deg = Right_deg; - deg_o = Left_deg; - new_memory(dim, deg, deg_o); - - Scaling_F = Scaling_F2 = 1.0; - - /*** y = 0 : ZY = LC ***/ - loop_ZY(dim, 0.0); - eval_frequency(dim, deg_o); - eval_Si_Si_1(dim, 0.0); - store_SiSv_1(dim, 0); - store(dim, 0); - - /*** Step 1 ***/ - /*** Step 2 ***/ - for (i = 1; i <= deg_o; i++) { - loop_ZY(dim, frequency[i]); - eval_Si_Si_1(dim, frequency[i]); - store_SiSv_1(dim, i); - store(dim, i); - } - poly_matrix(Sip, dim, deg_o); - poly_matrix(Si_1p, dim, deg_o); - poly_matrix(Sv_1p, dim, deg_o); - poly_W(dim, deg_o); - matrix_p_mult(Sip, W, Si_1p, dim, deg_o, deg_o, IWI); - matrix_p_mult(Sip, W, Sv_1p, dim, deg_o, deg_o, IWV); - - poly_matrix(SiSv_1, dim, deg_o); - - /*** - check(Sip, Si_1p, Sv_1p, SiSv_1); - ***/ - - generate_out(dim, deg_o); - - return(1); + int deg, deg_o; + int i; + + deg = Right_deg; + deg_o = Left_deg; + new_memory(dim, deg, deg_o); + + Scaling_F = Scaling_F2 = 1.0; + + /*** y = 0 : ZY = LC ***/ + loop_ZY(dim, 0.0); + eval_frequency(dim, deg_o); + eval_Si_Si_1(dim, 0.0); + store_SiSv_1(dim, 0); + store(dim, 0); + + /*** Step 1 ***/ + /*** Step 2 ***/ + for (i = 1; i <= deg_o; i++) { + loop_ZY(dim, frequency[i]); + eval_Si_Si_1(dim, frequency[i]); + store_SiSv_1(dim, i); + store(dim, i); + } + poly_matrix(Sip, dim, deg_o); + poly_matrix(Si_1p, dim, deg_o); + poly_matrix(Sv_1p, dim, deg_o); + poly_W(dim, deg_o); + matrix_p_mult(Sip, W, Si_1p, dim, deg_o, deg_o, IWI); + matrix_p_mult(Sip, W, Sv_1p, dim, deg_o, deg_o, IWV); + + poly_matrix(SiSv_1, dim, deg_o); + + /*** + check(Sip, Si_1p, Sv_1p, SiSv_1); + ***/ + + generate_out(dim, deg_o); + + return(1); } /*** @@ -1330,253 +1332,256 @@ coupled(int dim) static int generate_out(int dim, int deg_o) { - int i, j, k, rtv; - double C; - double *p; - double c1, c2, c3, x1, x2, x3; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - p = SiSv_1[i][j]; - SIV[i][j].C_0 = C = p[0]; - if (C == 0.0) - continue; - for (k = 0; k <= deg_o; k++) - p[k] /= C; - if (i == j) { - rtv = Pade_apx(sqrt(G_m[i][i] / R_m[i][i]) / C, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - SIV[i][j].C_0 = 0.0; - printf("SIV\n"); - continue; - } - } else { - rtv = Pade_apx(0.0, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - SIV[i][j].C_0 = 0.0; - printf("SIV\n"); - continue; - } - } - p = SIV[i][j].Poly = (double *) calloc(7, sizeof(double)); - p[0] = c1; - p[1] = c2; - p[2] = c3; - p[3] = x1; - p[4] = x2; - p[5] = x3; - p[6] = rtv; - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - for (k = 0; k < dim; k++) { - p = IWI[i][j].Poly[k]; - C = IWI[i][j].C_0[k]; - if (C == 0.0) - continue; - if (i == j && k == i) { - rtv = Pade_apx( - exp(- sqrt(G_m[i][i] * R_m[i][i]) * length) / C, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - IWI[i][j].C_0[k] = 0.0; - printf("IWI %d %d %d\n", i, j, k); - continue; - } - } else { - rtv = Pade_apx(0.0, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - IWI[i][j].C_0[k] = 0.0; - printf("IWI %d %d %d\n", i, j, k); - continue; - } - } - p[0] = c1; - p[1] = c2; - p[2] = c3; - p[3] = x1; - p[4] = x2; - p[5] = x3; - p[6] = rtv; - } - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - for (k = 0; k < dim; k++) { - p = IWV[i][j].Poly[k]; - C = IWV[i][j].C_0[k]; - if (C == 0.0) - continue; - if (i == j && k == i) { - rtv = Pade_apx(sqrt(G_m[i][i] / R_m[i][i]) * - exp(- sqrt(G_m[i][i] * R_m[i][i]) * length) / C, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - IWV[i][j].C_0[k] = 0.0; - printf("IWV %d %d %d\n", i, j, k); - continue; - } - } else { - rtv = Pade_apx(0.0, - p, &c1, &c2, &c3, &x1, &x2, &x3); - if (rtv == 0) { - IWV[i][j].C_0[k] = 0.0; - printf("IWV %d %d %d\n", i, j, k); - continue; - } - } - p[0] = c1; - p[1] = c2; - p[2] = c3; - p[3] = x1; - p[4] = x2; - p[5] = x3; - p[6] = rtv; - } - return(1); + int i, j, k, rtv; + double C; + double* p; + double c1, c2, c3, x1, x2, x3; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + p = SiSv_1[i][j]; + SIV[i][j].C_0 = C = p[0]; + if (C == 0.0) + continue; + for (k = 0; k <= deg_o; k++) + p[k] /= C; + if (i == j) { + rtv = Pade_apx(sqrt(G_m[i][i] / R_m[i][i]) / C, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + SIV[i][j].C_0 = 0.0; + printf("SIV\n"); + continue; + } + } + else { + rtv = Pade_apx(0.0, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + SIV[i][j].C_0 = 0.0; + printf("SIV\n"); + continue; + } + } + p = SIV[i][j].Poly = (double*)calloc(7, sizeof(double)); + p[0] = c1; + p[1] = c2; + p[2] = c3; + p[3] = x1; + p[4] = x2; + p[5] = x3; + p[6] = rtv; + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + for (k = 0; k < dim; k++) { + p = IWI[i][j].Poly[k]; + C = IWI[i][j].C_0[k]; + if (C == 0.0) + continue; + if (i == j && k == i) { + rtv = Pade_apx( + exp(-sqrt(G_m[i][i] * R_m[i][i]) * length) / C, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + IWI[i][j].C_0[k] = 0.0; + printf("IWI %d %d %d\n", i, j, k); + continue; + } + } + else { + rtv = Pade_apx(0.0, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + IWI[i][j].C_0[k] = 0.0; + printf("IWI %d %d %d\n", i, j, k); + continue; + } + } + p[0] = c1; + p[1] = c2; + p[2] = c3; + p[3] = x1; + p[4] = x2; + p[5] = x3; + p[6] = rtv; + } + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + for (k = 0; k < dim; k++) { + p = IWV[i][j].Poly[k]; + C = IWV[i][j].C_0[k]; + if (C == 0.0) + continue; + if (i == j && k == i) { + rtv = Pade_apx(sqrt(G_m[i][i] / R_m[i][i]) * + exp(-sqrt(G_m[i][i] * R_m[i][i]) * length) / C, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + IWV[i][j].C_0[k] = 0.0; + printf("IWV %d %d %d\n", i, j, k); + continue; + } + } + else { + rtv = Pade_apx(0.0, + p, &c1, &c2, &c3, &x1, &x2, &x3); + if (rtv == 0) { + IWV[i][j].C_0[k] = 0.0; + printf("IWV %d %d %d\n", i, j, k); + continue; + } + } + p[0] = c1; + p[1] = c2; + p[2] = c3; + p[3] = x1; + p[4] = x2; + p[5] = x3; + p[6] = rtv; + } + return(1); } /**************************************************************** - mult.c Multiplication for Matrix of Polynomial - X(y) = A(y) D(y) B(y), - where D(y) is a diagonal matrix with each - diagonal entry of the form - e^{-a_i s}d(y), for which s = 1/y - and i = 1..N. - Each entry of X(y) will be of the form - \sum_{i=1}^N c_i e^{-a_i s} b_i(y), where - b_i(0) = 1; therefore, those - b_i(y)'s will be each entry's output. + mult.c Multiplication for Matrix of Polynomial + X(y) = A(y) D(y) B(y), + where D(y) is a diagonal matrix with each + diagonal entry of the form + e^{-a_i s}d(y), for which s = 1/y + and i = 1..N. + Each entry of X(y) will be of the form + \sum_{i=1}^N c_i e^{-a_i s} b_i(y), where + b_i(0) = 1; therefore, those + b_i(y)'s will be each entry's output. ****************************************************************/ static void -mult_p(double *p1, double *p2, double *p3, int d1, int d2, int d3) +mult_p(double* p1, double* p2, double* p3, int d1, int d2, int d3) /* p3 = p1 * p2 */ { - int i, j, k; - - for (i = 0; i <= d3; i++) - p3[i] = 0.0; - for (i = 0; i <= d1; i++) - for (j = i, k = 0; k <= d2; j++, k++) { - if (j > d3) - break; - p3[j] += p1[i] * p2[k]; - } + int i, j, k; + + for (i = 0; i <= d3; i++) + p3[i] = 0.0; + for (i = 0; i <= d1; i++) + for (j = i, k = 0; k <= d2; j++, k++) { + if (j > d3) + break; + p3[j] += p1[i] * p2[k]; + } } static void matrix_p_mult( - double *A_in[MAX_DIM][MAX_DIM], - double *D1[MAX_DIM], - double *B[MAX_DIM][MAX_DIM], - int dim, int deg, int deg_o, - Mult_Out X[MAX_DIM][MAX_DIM]) + double* A_in[MAX_DIM][MAX_DIM], + double* D1[MAX_DIM], + double* B[MAX_DIM][MAX_DIM], + int dim, int deg, int deg_o, + Mult_Out X[MAX_DIM][MAX_DIM]) { - int i, j, k, l; - double *p; - double *T[MAX_DIM][MAX_DIM]; - double t1; - - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - p = T[i][j] = (double *) calloc((size_t) (deg_o+1), sizeof(double)); - mult_p(B[i][j], D1[i], p, deg, deg_o, deg_o); - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - for (k = 0; k < dim; k++) { - p = X[i][j].Poly[k] = - (double *) calloc((size_t) (deg_o+1), sizeof(double)); - mult_p(A_in[i][k], T[k][j], p, deg, deg_o, deg_o); - t1 = X[i][j].C_0[k] = p[0]; - if (t1 != 0.0) { - p[0] = 1.0; - for (l = 1; l <= deg_o; l++) - p[l] /= t1; - } - } - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) - tfree(T[i][j]); - - /********** - for (i = 0; i < dim; i++) - for (j = 0; j < dim; j++) { - for (k = 0; k < dim; k++) { - fprintf(outFile, "(%5.3f)", X[i][j].C_0[k]); - p = X[i][j].Poly[k]; - for (l = 0; l <= deg_o; l++) - fprintf(outFile, "%5.3f ", p[l]); - fprintf(outFile, "\n"); - } - fprintf(outFile, "\n"); - } - ***********/ + int i, j, k, l; + double* p; + double* T[MAX_DIM][MAX_DIM]; + double t1; + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + p = T[i][j] = (double*)calloc((size_t)(deg_o + 1), sizeof(double)); + mult_p(B[i][j], D1[i], p, deg, deg_o, deg_o); + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + for (k = 0; k < dim; k++) { + p = X[i][j].Poly[k] = + (double*)calloc((size_t)(deg_o + 1), sizeof(double)); + mult_p(A_in[i][k], T[k][j], p, deg, deg_o, deg_o); + t1 = X[i][j].C_0[k] = p[0]; + if (t1 != 0.0) { + p[0] = 1.0; + for (l = 1; l <= deg_o; l++) + p[l] /= t1; + } + } + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + tfree(T[i][j]); + + /********** + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) { + for (k = 0; k < dim; k++) { + fprintf(outFile, "(%5.3f)", X[i][j].C_0[k]); + p = X[i][j].Poly[k]; + for (l = 0; l <= deg_o; l++) + fprintf(outFile, "%5.3f ", p[l]); + fprintf(outFile, "\n"); + } + fprintf(outFile, "\n"); + } + ***********/ } /**************************************************************** - mode approximation + mode approximation ****************************************************************/ -/*** - ***/ + /*** + ***/ static double -approx_mode(double *X, double *b, double length_in) +approx_mode(double* X, double* b, double length_in) { - double w0, w1, w2, w3, w4, w5; - double a[8]; - double delay, atten; - double y1, y2, y3, y4, y5, y6; - int i, j; - - w0 = X[0]; - w1 = X[1] / w0; /* a */ - w2 = X[2] / w0; /* b */ - w3 = X[3] / w0; /* c */ - w4 = X[4] / w0; /* d */ - w5 = X[5] / w0; /* e */ - - y1 = 0.5 * w1; - y2 = w2 - y1 * y1; - y3 = 3 * w3 - 3.0 * y1 * y2; - y4 = 12.0 * w4 - 3.0 * y2 * y2 - 4.0 * y1 * y3; - y5 = 60.0 * w5 - 5.0 * y1 * y4 -10.0 * y2 * y3; - y6 = -10.0 * y3 * y3 - 15.0 * y2 * y4 - 6.0 * y1 * y5; - - delay = sqrt(w0) * length_in / Scaling_F; - atten = exp(- delay * y1); - - a[1] = y2 / 2.0; - a[2] = y3 / 6.0; - a[3] = y4 / 24.0; - a[4] = y5 / 120.0; - a[5] = y6 / 720.0; - - a[1] *= -delay; - a[2] *= -delay; - a[3] *= -delay; - a[4] *= -delay; - a[5] *= -delay; - - b[0] = 1.0; - b[1] = a[1]; - for (i = 2; i <= 5; i++) { - b[i] = 0.0; - for (j = 1; j <= i; j++) - b[i] += j * a[j] * b[i-j]; - b[i] = b[i] / (double) i; - } - - for (i = 0; i <= 5; i++) - b[i] *= atten; - - return(delay); + double w0, w1, w2, w3, w4, w5; + double a[8]; + double delay, atten; + double y1, y2, y3, y4, y5, y6; + int i, j; + + w0 = X[0]; + w1 = X[1] / w0; /* a */ + w2 = X[2] / w0; /* b */ + w3 = X[3] / w0; /* c */ + w4 = X[4] / w0; /* d */ + w5 = X[5] / w0; /* e */ + + y1 = 0.5 * w1; + y2 = w2 - y1 * y1; + y3 = 3 * w3 - 3.0 * y1 * y2; + y4 = 12.0 * w4 - 3.0 * y2 * y2 - 4.0 * y1 * y3; + y5 = 60.0 * w5 - 5.0 * y1 * y4 - 10.0 * y2 * y3; + y6 = -10.0 * y3 * y3 - 15.0 * y2 * y4 - 6.0 * y1 * y5; + + delay = sqrt(w0) * length_in / Scaling_F; + atten = exp(-delay * y1); + + a[1] = y2 / 2.0; + a[2] = y3 / 6.0; + a[3] = y4 / 24.0; + a[4] = y5 / 120.0; + a[5] = y6 / 720.0; + + a[1] *= -delay; + a[2] *= -delay; + a[3] *= -delay; + a[4] *= -delay; + a[5] *= -delay; + + b[0] = 1.0; + b[1] = a[1]; + for (i = 2; i <= 5; i++) { + b[i] = 0.0; + for (j = 1; j <= i; j++) + b[i] += j * a[j] * b[i - j]; + b[i] = b[i] / (double)i; + } + + for (i = 0; i <= 5; i++) + b[i] *= atten; + + return(delay); } /*** @@ -1585,7 +1590,7 @@ approx_mode(double *X, double *b, double length_in) static double eval2(double a, double b, double c, double x) { - return(a*x*x + b*x + c); + return(a * x * x + b * x + c); } /*** @@ -1593,347 +1598,354 @@ eval2(double a, double b, double c, double x) static int get_c(double q1, double q2, double q3, double p1, double p2, double a, double b, - double *cr, double *ci) + double* cr, double* ci) { - double d, n; - - d = (3.0*(a*a-b*b)+2.0*p1*a+p2)*(3.0*(a*a-b*b)+2.0*p1*a+p2); - d += (6.0*a*b+2.0*p1*b)*(6.0*a*b+2.0*p1*b); - n = -(q1*(a*a-b*b)+q2*a+q3)*(6.0*a*b+2.0*p1*b); - n += (2.0*q1*a*b+q2*b)*(3.0*(a*a-b*b)+2.0*p1*a+p2); - *ci = n/d; - n = (3.0*(a*a-b*b)+2.0*p1*a+p2)*(q1*(a*a-b*b)+q2*a+q3); - n += (6.0*a*b+2.0*p1*b)*(2.0*q1*a*b+q2*b); - *cr = n/d; - - return(1); + double d, n; + + d = (3.0 * (a * a - b * b) + 2.0 * p1 * a + p2) * (3.0 * (a * a - b * b) + 2.0 * p1 * a + p2); + d += (6.0 * a * b + 2.0 * p1 * b) * (6.0 * a * b + 2.0 * p1 * b); + n = -(q1 * (a * a - b * b) + q2 * a + q3) * (6.0 * a * b + 2.0 * p1 * b); + n += (2.0 * q1 * a * b + q2 * b) * (3.0 * (a * a - b * b) + 2.0 * p1 * a + p2); + *ci = n / d; + n = (3.0 * (a * a - b * b) + 2.0 * p1 * a + p2) * (q1 * (a * a - b * b) + q2 * a + q3); + n += (6.0 * a * b + 2.0 * p1 * b) * (2.0 * q1 * a * b + q2 * b); + *cr = n / d; + + return(1); } static int -Pade_apx(double a_b, double *b, double *c1, double *c2, double *c3, - double *x1, double *x2, double *x3) -/* - b[0] + b[1]*y + b[2]*y^2 + ... + b[5]*y^5 + ... - = (q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1) - - where b[0] is always equal to 1.0 and neglected, - and y = 1/s. - - (q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1) - = (s^3 + q1*s^2 + q2*s + q3) / (s^3 + p1*s^2 + p2*s + p3) - = c1 / (s - x1) + c2 / (s - x2) + c3 / (s - x3) + 1.0 - */ +Pade_apx(double a_b, double* b, double* c1, double* c2, double* c3, + double* x1, double* x2, double* x3) + /* + b[0] + b[1]*y + b[2]*y^2 + ... + b[5]*y^5 + ... + = (q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1) + + where b[0] is always equal to 1.0 and neglected, + and y = 1/s. + + (q3*y^3 + q2*y^2 + q1*y + 1) / (p3*y^3 + p2*y^2 + p1*y + 1) + = (s^3 + q1*s^2 + q2*s + q3) / (s^3 + p1*s^2 + p2*s + p3) + = c1 / (s - x1) + c2 / (s - x2) + c3 / (s - x3) + 1.0 + */ { - double p1, p2, p3, q1, q2, q3; - - At[0][0] = 1.0 - a_b; - At[0][1] = b[1]; - At[0][2] = b[2]; - At[0][3] = -b[3]; - - At[1][0] = b[1]; - At[1][1] = b[2]; - At[1][2] = b[3]; - At[1][3] = -b[4]; - - At[2][0] = b[2]; - At[2][1] = b[3]; - At[2][2] = b[4]; - At[2][3] = -b[5]; - - Gaussian_Elimination(3); - - p3 = At[0][3]; - p2 = At[1][3]; - p1 = At[2][3]; - /* - if (p3 < 0.0 || p2 < 0.0 || p1 < 0.0 || p1*p2 <= p3) - return(0); - */ - q1 = p1 + b[1]; - q2 = b[1] * p1 + p2 + b[2]; - q3 = p3 * a_b; - - if (find_roots(p1, p2, p3, x1, x2, x3)) { - /* - printf("complex roots : %e %e %e \n", *x1, *x2, *x3); - */ - *c1 = eval2(q1 - p1, q2 - p2, q3 - p3, *x1) / - eval2(3.0, 2.0 * p1, p2, *x1); - get_c(q1 - p1, q2 - p2, q3 - p3, p1, p2, *x2, *x3, c2, c3); - return(2); - } else { - /* new - printf("roots are %e %e %e \n", *x1, *x2, *x3); - */ - *c1 = eval2(q1 - p1, q2 - p2, q3 - p3, *x1) / - eval2(3.0, 2.0 * p1, p2, *x1); - *c2 = eval2(q1 - p1, q2 - p2, q3 - p3, *x2) / - eval2(3.0, 2.0 * p1, p2, *x2); - *c3 = eval2(q1 - p1, q2 - p2, q3 - p3, *x3) / - eval2(3.0, 2.0 * p1, p2, *x3); - return(1); - } + double p1, p2, p3, q1, q2, q3; + + At[0][0] = 1.0 - a_b; + At[0][1] = b[1]; + At[0][2] = b[2]; + At[0][3] = -b[3]; + + At[1][0] = b[1]; + At[1][1] = b[2]; + At[1][2] = b[3]; + At[1][3] = -b[4]; + + At[2][0] = b[2]; + At[2][1] = b[3]; + At[2][2] = b[4]; + At[2][3] = -b[5]; + + Gaussian_Elimination(3); + + p3 = At[0][3]; + p2 = At[1][3]; + p1 = At[2][3]; + /* + if (p3 < 0.0 || p2 < 0.0 || p1 < 0.0 || p1*p2 <= p3) + return(0); + */ + q1 = p1 + b[1]; + q2 = b[1] * p1 + p2 + b[2]; + q3 = p3 * a_b; + + if (find_roots(p1, p2, p3, x1, x2, x3)) { + /* + printf("complex roots : %e %e %e \n", *x1, *x2, *x3); + */ + *c1 = eval2(q1 - p1, q2 - p2, q3 - p3, *x1) / + eval2(3.0, 2.0 * p1, p2, *x1); + get_c(q1 - p1, q2 - p2, q3 - p3, p1, p2, *x2, *x3, c2, c3); + return(2); + } + else { + /* new + printf("roots are %e %e %e \n", *x1, *x2, *x3); + */ + *c1 = eval2(q1 - p1, q2 - p2, q3 - p3, *x1) / + eval2(3.0, 2.0 * p1, p2, *x1); + *c2 = eval2(q1 - p1, q2 - p2, q3 - p3, *x2) / + eval2(3.0, 2.0 * p1, p2, *x2); + *c3 = eval2(q1 - p1, q2 - p2, q3 - p3, *x3) / + eval2(3.0, 2.0 * p1, p2, *x3); + return(1); + } } static int Gaussian_Elimination(int dims) { - int i, j, k, dim; - double f; - double max; - int imax; - - dim = dims; - - for (i = 0; i < dim; i++) { - imax = i; - max = ABS(At[i][i]); - for (j = i+1; j < dim; j++) - if (ABS(At[j][i]) > max) { - imax = j; - max = ABS(At[j][i]); - } - if (max < epsi_mult) { - fprintf(stderr, " can not choose a pivot (mult)\n"); - controlled_exit(EXIT_FAILURE); - } - if (imax != i) - for (k = i; k <= dim; k++) { - SWAP(double, At[i][k], At[imax][k]); - } - - f = 1.0 / At[i][i]; - At[i][i] = 1.0; - - for (j = i+1; j <= dim; j++) - At[i][j] *= f; - - for (j = 0; j < dim ; j++) { - if (i == j) - continue; - f = At[j][i]; - At[j][i] = 0.0; - for (k = i+1; k <= dim; k++) - At[j][k] -= f * At[i][k]; - } - } - return(1); + int i, j, k, dim; + double f; + double max; + int imax; + + dim = dims; + + for (i = 0; i < dim; i++) { + imax = i; + max = ABS(At[i][i]); + for (j = i + 1; j < dim; j++) + if (ABS(At[j][i]) > max) { + imax = j; + max = ABS(At[j][i]); + } + if (max < epsi_mult) { + fprintf(stderr, " can not choose a pivot (mult)\n"); + controlled_exit(EXIT_FAILURE); + } + if (imax != i) + for (k = i; k <= dim; k++) { + SWAP(double, At[i][k], At[imax][k]); + } + + f = 1.0 / At[i][i]; + At[i][i] = 1.0; + + for (j = i + 1; j <= dim; j++) + At[i][j] *= f; + + for (j = 0; j < dim; j++) { + if (i == j) + continue; + f = At[j][i]; + At[j][i] = 0.0; + for (k = i + 1; k <= dim; k++) + At[j][k] -= f * At[i][k]; + } + } + return(1); } static double root3(double a1, double a2, double a3, double x) { - double t1, t2; + double t1, t2; - t1 = x * (x * (x + a1) + a2) + a3; - t2 = x * (2.0*a1 + 3.0*x) + a2; + t1 = x * (x * (x + a1) + a2) + a3; + t2 = x * (2.0 * a1 + 3.0 * x) + a2; - return(x - t1 / t2); + return(x - t1 / t2); } static int -div3(double a1, double a2, double a3, double x, double *p1, double *p2) +div3(double a1, double a2, double a3, double x, double* p1, double* p2) { - NG_IGNORE(a2); + NG_IGNORE(a2); - *p1 = a1 + x; + *p1 = a1 + x; - /* *p2 = a2 + (a1 + x) * x; */ + /* *p2 = a2 + (a1 + x) * x; */ - *p2 = - a3 / x; + *p2 = -a3 / x; - return(1); + return(1); } static int -find_roots(double a1, double a2, double a3, double *x1, double *x2, double *x3) +find_roots(double a1, double a2, double a3, double* x1, double* x2, double* x3) { - double x, t; - double p, q; - - /*********************************************** - double m,n; - p = a1*a1/3.0 - a2; - q = a1*a2/3.0 - a3 - 2.0*a1*a1*a1/27.0; - p = p*p*p/27.0; - t = q*q - 4.0*p; - if (t < 0.0) { - if (q != 0.0) { - t = atan(sqrt((double)-t)/q); - if (t < 0.0) - t += M_PI; - t /= 3.0; - x = 2.0 * pow(p, 0.16666667) * cos(t) - a1 / 3.0; - } else { - t /= -4.0; - x = pow(t, 0.16666667) * 1.732 - a1 / 3.0; - } - } else { - t = sqrt(t); - m = 0.5*(q - t); - n = 0.5*(q + t); - if (m < 0.0) - m = -pow((double) -m, (double) 0.3333333); - else - m = pow((double) m, (double) 0.3333333); - if (n < 0.0) - n = -pow((double) -n, (double) 0.3333333); - else - n = pow((double) n, (double) 0.3333333); - x = m + n - a1 / 3.0; - } - ************************************************/ - q = (a1*a1-3.0*a2) / 9.0; - p = (2.0*a1*a1*a1-9.0*a1*a2+27.0*a3) / 54.0; - t = q*q*q - p*p; - if (t >= 0.0) { - t = acos(p /(q * sqrt(q))); - x = -2.0*sqrt(q)*cos(t / 3.0) - a1/3.0; - } else { - if (p > 0.0) { - t = pow(sqrt(-t)+p, 1.0 / 3.0); - x = -(t + q / t) - a1/3.0; - } else if (p == 0.0) { - x = -a1/3.0; - } else { - t = pow(sqrt(-t)-p, 1.0 / 3.0); - x = (t + q / t) - a1/3.0; - } - } - /* - fprintf(stderr, "..1.. %e\n", x*x*x+a1*x*x+a2*x+a3); - */ - { - double x_backup = x; - int i = 0; - for (t = root3(a1, a2, a3, x); ABS(t-x) > 5.0e-4; - t = root3(a1, a2, a3, x)) - if (++i == 32) { - x = x_backup; - break; - } else - x = t; - } - /* - fprintf(stderr, "..2.. %e\n", x*x*x+a1*x*x+a2*x+a3); - */ - - - *x1 = x; - div3(a1, a2, a3, x, &a1, &a2); - - t = a1 * a1 - 4.0 * a2; - if (t < 0) { - /* - fprintf(stderr, "***** Two Imaginary Roots.\n Update.\n"); - *x2 = -0.5 * a1; - *x3 = a2 / *x2; - */ - *x3 = 0.5 * sqrt(-t); - *x2 = -0.5 * a1; - return(1); - } else { - t = sqrt(t); - if (a1 >= 0.0) - *x2 = t = -0.5 * (a1 + t); - else - *x2 = t = -0.5 * (a1 - t); - *x3 = a2 / t; - return(0); - } + double x, t; + double p, q; + + /*********************************************** + double m,n; + p = a1*a1/3.0 - a2; + q = a1*a2/3.0 - a3 - 2.0*a1*a1*a1/27.0; + p = p*p*p/27.0; + t = q*q - 4.0*p; + if (t < 0.0) { + if (q != 0.0) { + t = atan(sqrt((double)-t)/q); + if (t < 0.0) + t += M_PI; + t /= 3.0; + x = 2.0 * pow(p, 0.16666667) * cos(t) - a1 / 3.0; + } else { + t /= -4.0; + x = pow(t, 0.16666667) * 1.732 - a1 / 3.0; + } + } else { + t = sqrt(t); + m = 0.5*(q - t); + n = 0.5*(q + t); + if (m < 0.0) + m = -pow((double) -m, (double) 0.3333333); + else + m = pow((double) m, (double) 0.3333333); + if (n < 0.0) + n = -pow((double) -n, (double) 0.3333333); + else + n = pow((double) n, (double) 0.3333333); + x = m + n - a1 / 3.0; + } + ************************************************/ + q = (a1 * a1 - 3.0 * a2) / 9.0; + p = (2.0 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0; + t = q * q * q - p * p; + if (t >= 0.0) { + t = acos(p / (q * sqrt(q))); + x = -2.0 * sqrt(q) * cos(t / 3.0) - a1 / 3.0; + } + else { + if (p > 0.0) { + t = pow(sqrt(-t) + p, 1.0 / 3.0); + x = -(t + q / t) - a1 / 3.0; + } + else if (p == 0.0) { + x = -a1 / 3.0; + } + else { + t = pow(sqrt(-t) - p, 1.0 / 3.0); + x = (t + q / t) - a1 / 3.0; + } + } + /* + fprintf(stderr, "..1.. %e\n", x*x*x+a1*x*x+a2*x+a3); + */ + { + double x_backup = x; + int i = 0; + for (t = root3(a1, a2, a3, x); ABS(t - x) > 5.0e-4; + t = root3(a1, a2, a3, x)) + if (++i == 32) { + x = x_backup; + break; + } + else + x = t; + } + /* + fprintf(stderr, "..2.. %e\n", x*x*x+a1*x*x+a2*x+a3); + */ + + + *x1 = x; + div3(a1, a2, a3, x, &a1, &a2); + + t = a1 * a1 - 4.0 * a2; + if (t < 0) { + /* + fprintf(stderr, "***** Two Imaginary Roots.\n Update.\n"); + *x2 = -0.5 * a1; + *x3 = a2 / *x2; + */ + *x3 = 0.5 * sqrt(-t); + *x2 = -0.5 * a1; + return(1); + } + else { + t = sqrt(t); + if (a1 >= 0.0) + *x2 = t = -0.5 * (a1 + t); + else + *x2 = t = -0.5 * (a1 - t); + *x3 = a2 / t; + return(0); + } } static NDnamePt -insert_ND(char *name, NDnamePt *ndn) +insert_ND(char* name, NDnamePt* ndn) { - int cmp; - NDnamePt p; - - if (*ndn == NULL) { - p = *ndn = TMALLOC(NDname, 1); - p->nd = NULL; - p->right = p->left = NULL; - strcpy(p->id, name); - return(p); - } - cmp = strcmp((*ndn)->id, name); - if (cmp == 0) - return(*ndn); - else { - if (cmp < 0) - return(insert_ND(name, &((*ndn)->left))); - else - return(insert_ND(name, &((*ndn)->right))); - } + int cmp; + NDnamePt p; + + if (*ndn == NULL) { + p = *ndn = TMALLOC(NDname, 1); + p->nd = NULL; + p->right = p->left = NULL; + strcpy(p->id, name); + return(p); + } + cmp = strcmp((*ndn)->id, name); + if (cmp == 0) + return(*ndn); + else { + if (cmp < 0) + return(insert_ND(name, &((*ndn)->left))); + else + return(insert_ND(name, &((*ndn)->right))); + } } -static NODE * -insert_node(char *name) +static NODE* +insert_node(char* name) { - NDnamePt n; - NODE *p; - - n = insert_ND(name, &ndn_btree); - if (n->nd == NULL) { - p = NEW_node(); - p->name = n; - n->nd = p; - p->next = node_tab; - node_tab = p; - return(p); - } else - return(n->nd); + NDnamePt n; + NODE* p; + + n = insert_ND(name, &ndn_btree); + if (n->nd == NULL) { + p = NEW_node(); + p->name = n; + n->nd = p; + p->next = node_tab; + node_tab = p; + return(p); + } + else + return(n->nd); } /*** static int divC(double ar, double ai, double br, double bi, double *cr, double *ci) { - double t; + double t; - t = br*br + bi*bi; - *cr = (ar*br + ai*bi) / t; - *ci = (ai*br - ar*bi) / t; + t = br*br + bi*bi; + *cr = (ar*br + ai*bi) / t; + *ci = (ai*br - ar*bi) / t; - return(1); + return(1); } ***/ static NODE -*NEW_node(void) +* NEW_node(void) { - NODE *n; - - n = TMALLOC(NODE, 1); - n->mptr = NULL; - n->gptr = NULL; - n->cptr = NULL; - n->rptr = NULL; - n->tptr = NULL; - n->cplptr = NULL; - n->rlptr = NULL; - n->ddptr = NULL; - n->cvccsptr = NULL; - n->vccsptr = NULL; - n->CL = 0.001; - n->V = n->dv = 0.0; - n->gsum = n->cgsum = 0; - n->is = 0; - n->tag = 0; - n->flag = 0; - n->region = NULL; - n->ofile = NULL; - n->dvtag = 0; - - return(n); + NODE* n; + + n = TMALLOC(NODE, 1); + n->mptr = NULL; + n->gptr = NULL; + n->cptr = NULL; + n->rptr = NULL; + n->tptr = NULL; + n->cplptr = NULL; + n->rlptr = NULL; + n->ddptr = NULL; + n->cvccsptr = NULL; + n->vccsptr = NULL; + n->CL = 0.001; + n->V = n->dv = 0.0; + n->gsum = n->cgsum = 0; + n->is = 0; + n->tag = 0; + n->flag = 0; + n->region = NULL; + n->ofile = NULL; + n->dvtag = 0; + + return(n); } /**************************************************************** - diag.c This file contains the main(). + diag.c This file contains the main(). ****************************************************************/ #define epsi2 1.0e-8 @@ -1944,210 +1956,211 @@ static MAXE_PTR row; static MAXE_PTR sort(MAXE_PTR list, double val, int r, int c, MAXE_PTR e) { - if (list == NULL || list->value < val) { - e->row = r; - e->col = c; - e->value = val; - e->next = list; - return(e); - } else { - list->next = sort(list->next, val, r, c, e); - return(list); - } + if (list == NULL || list->value < val) { + e->row = r; + e->col = c; + e->value = val; + e->next = list; + return(e); + } + else { + list->next = sort(list->next, val, r, c, e); + return(list); + } } static void ordering(void) { - MAXE_PTR e; - int i, j, m; - double mv; - - for (i = 0; i < dim-1; i++) { - m = i+1; - mv = ABS(ZY[i][m]); - for (j = m+1; j < dim; j++) - if ((int)(ABS(ZY[i][j]) * 1e7) > (int) (1e7 *mv)) { - - mv = ABS(ZY[i][j]); - m = j; - } - e = TMALLOC(MAXE, 1); - row = sort(row, mv, i, m, e); - } + MAXE_PTR e; + int i, j, m; + double mv; + + for (i = 0; i < dim - 1; i++) { + m = i + 1; + mv = ABS(ZY[i][m]); + for (j = m + 1; j < dim; j++) + if ((int)(ABS(ZY[i][j]) * 1e7) > (int) (1e7 * mv)) { + + mv = ABS(ZY[i][j]); + m = j; + } + e = TMALLOC(MAXE, 1); + row = sort(row, mv, i, m, e); + } } static MAXE_PTR -delete_1(MAXE_PTR *list, int rc) +delete_1(MAXE_PTR* list, int rc) { - MAXE_PTR list1, e; - - list1 = *list; - if ((*list)->row == rc) { - *list = (*list)->next; - return(list1); - } - for (e = list1->next; e->row != rc; e = e->next) - list1 = e; - list1->next = e->next; - return(e); + MAXE_PTR list1, e; + + list1 = *list; + if ((*list)->row == rc) { + *list = (*list)->next; + return(list1); + } + for (e = list1->next; e->row != rc; e = e->next) + list1 = e; + list1->next = e->next; + return(e); } static void reordering(int p, int q) { - MAXE_PTR e; - int j, m; - double mv; - - m = p+1; - mv = ABS(ZY[p][m]); - for (j = m+1; j < dim; j++) - if ((int)(ABS(ZY[p][j]) * 1e7) > (int) (1e7 *mv)) { - mv = ABS(ZY[p][j]); - m = j; - } - e = delete_1(&row, p); - row = sort(row, mv, p, m, e); - - m = q+1; - if (m != dim) { - mv = ABS(ZY[q][m]); - for (j = m+1; j < dim; j++) - if ((int)(ABS(ZY[q][j]) * 1e7) > (int) (1e7 *mv)) { - - mv = ABS(ZY[q][j]); - m = j; - } - e = delete_1(&row, q); - row = sort(row, mv, q, m, e); - } + MAXE_PTR e; + int j, m; + double mv; + + m = p + 1; + mv = ABS(ZY[p][m]); + for (j = m + 1; j < dim; j++) + if ((int)(ABS(ZY[p][j]) * 1e7) > (int) (1e7 * mv)) { + mv = ABS(ZY[p][j]); + m = j; + } + e = delete_1(&row, p); + row = sort(row, mv, p, m, e); + + m = q + 1; + if (m != dim) { + mv = ABS(ZY[q][m]); + for (j = m + 1; j < dim; j++) + if ((int)(ABS(ZY[q][j]) * 1e7) > (int) (1e7 * mv)) { + + mv = ABS(ZY[q][j]); + m = j; + } + e = delete_1(&row, q); + row = sort(row, mv, q, m, e); + } } static void diag(int dims) { - int i, j, c; - double fmin, fmax; - - dim = dims; - row = NULL; - - fmin = fmax = ABS(ZY[0][0]); - for (i = 0; i < dim; i++) - for (j = i; j < dim; j++) - if (ABS(ZY[i][j]) > fmax) - fmax = ABS(ZY[i][j]); - else if (ABS(ZY[i][j]) < fmin) - fmin = ABS(ZY[i][j]); - fmin = 2.0 / (fmin + fmax); - for (i = 0; i < dim; i++) - for (j = i; j < dim; j++) - ZY[i][j] *= fmin; - - for (i = 0; i < dim; i++) { - for (j = 0; j < dim; j++) - if (i == j) - Sv[i][i] = 1.0; - else - Sv[i][j] = 0.0; - } - - ordering(); - - if (row) - for (c = 0; row->value > epsi2; c++) { - int p, q; - - p = row->row; - q = row->col; - - rotate(dim, p, q); - reordering(p, q); - } - - for (i = 0; i < dim; i++) - D[i] = ZY[i][i] / fmin; - - while (row) { - MAXE_PTR tmp_row = row->next; - tfree(row); - row = tmp_row; - } + int i, j, c; + double fmin, fmax; + + dim = dims; + row = NULL; + + fmin = fmax = ABS(ZY[0][0]); + for (i = 0; i < dim; i++) + for (j = i; j < dim; j++) + if (ABS(ZY[i][j]) > fmax) + fmax = ABS(ZY[i][j]); + else if (ABS(ZY[i][j]) < fmin) + fmin = ABS(ZY[i][j]); + fmin = 2.0 / (fmin + fmax); + for (i = 0; i < dim; i++) + for (j = i; j < dim; j++) + ZY[i][j] *= fmin; + + for (i = 0; i < dim; i++) { + for (j = 0; j < dim; j++) + if (i == j) + Sv[i][i] = 1.0; + else + Sv[i][j] = 0.0; + } + + ordering(); + + if (row) + for (c = 0; row->value > epsi2; c++) { + int p, q; + + p = row->row; + q = row->col; + + rotate(dim, p, q); + reordering(p, q); + } + + for (i = 0; i < dim; i++) + D[i] = ZY[i][i] / fmin; + + while (row) { + MAXE_PTR tmp_row = row->next; + tfree(row); + row = tmp_row; + } } /**************************************************************** - rotate() rotation of the Jacobi's method + rotate() rotation of the Jacobi's method ****************************************************************/ static int rotate(int dim_in, int p, int q) { - int j; - double co, si; - double ve, mu, ld; - double T[MAX_DIM]; - double t; - - ld = - ZY[p][q]; - mu = 0.5 * (ZY[p][p] - ZY[q][q]); - ve = sqrt(ld*ld + mu*mu); - co = sqrt((ve + ABS(mu)) / (2.0 * ve)); - si = SGN(mu) * ld / (2.0 * ve * co); - - for (j = p+1; j < dim_in; j++) - T[j] = ZY[p][j]; - for (j = 0; j < p; j++) - T[j] = ZY[j][p]; - - for (j = p+1; j < dim_in; j++) { - if (j == q) - continue; - if (j > q) - ZY[p][j] = T[j] * co - ZY[q][j] * si; - else - ZY[p][j] = T[j] * co - ZY[j][q] * si; - } - for (j = q+1; j < dim_in; j++) { - if (j == p) - continue; - ZY[q][j] = T[j] * si + ZY[q][j] * co; - } - for (j = 0; j < p; j++) { - if (j == q) - continue; - ZY[j][p] = T[j] * co - ZY[j][q] * si; - } - for (j = 0; j < q; j++) { - if (j == p) - continue; - ZY[j][q] = T[j] * si + ZY[j][q] * co; - } - - t = ZY[p][p]; - ZY[p][p] = t * co * co + ZY[q][q] * si * si - 2.0 * ZY[p][q] * si * co; - ZY[q][q] = t * si * si + ZY[q][q] * co * co + 2.0 * ZY[p][q] * si * co; - - ZY[p][q] = 0.0; - - { - double R[MAX_DIM]; - - for (j = 0; j < dim_in; j++) { - T[j] = Sv[j][p]; - R[j] = Sv[j][q]; - } - - for (j = 0; j < dim_in; j++) { - Sv[j][p] = T[j] * co - R[j] * si; - Sv[j][q] = T[j] * si + R[j] * co; - } - } - - return(1); + int j; + double co, si; + double ve, mu, ld; + double T[MAX_DIM]; + double t; + + ld = -ZY[p][q]; + mu = 0.5 * (ZY[p][p] - ZY[q][q]); + ve = sqrt(ld * ld + mu * mu); + co = sqrt((ve + ABS(mu)) / (2.0 * ve)); + si = SGN(mu) * ld / (2.0 * ve * co); + + for (j = p + 1; j < dim_in; j++) + T[j] = ZY[p][j]; + for (j = 0; j < p; j++) + T[j] = ZY[j][p]; + + for (j = p + 1; j < dim_in; j++) { + if (j == q) + continue; + if (j > q) + ZY[p][j] = T[j] * co - ZY[q][j] * si; + else + ZY[p][j] = T[j] * co - ZY[j][q] * si; + } + for (j = q + 1; j < dim_in; j++) { + if (j == p) + continue; + ZY[q][j] = T[j] * si + ZY[q][j] * co; + } + for (j = 0; j < p; j++) { + if (j == q) + continue; + ZY[j][p] = T[j] * co - ZY[j][q] * si; + } + for (j = 0; j < q; j++) { + if (j == p) + continue; + ZY[j][q] = T[j] * si + ZY[j][q] * co; + } + + t = ZY[p][p]; + ZY[p][p] = t * co * co + ZY[q][q] * si * si - 2.0 * ZY[p][q] * si * co; + ZY[q][q] = t * si * si + ZY[q][q] * co * co + 2.0 * ZY[p][q] * si * co; + + ZY[p][q] = 0.0; + + { + double R[MAX_DIM]; + + for (j = 0; j < dim_in; j++) { + T[j] = Sv[j][p]; + R[j] = Sv[j][q]; + } + + for (j = 0; j < dim_in; j++) { + Sv[j][p] = T[j] * co - R[j] * si; + Sv[j][q] = T[j] * si + R[j] * co; + } + } + + return(1); }