/*********************************************************************
 *  FUNCTION:  calc_res -- We're trying to solve    MU = K           *
 *                         This function calculates the residual r:  *
 *                                                                   *
 *                                 R = K - MU                        *
 *                                                                   *
 *  INPUTS: R, M, U, K                                               *
 *********************************************************************/
void calc_res (R, M, U, K)
struct collvec *R, *U, *K;
struct collmat *M;
{
struct kvec2 tempk1, tempk2;
struct kvec4 k41, k42;
int i, im1;

/***************************************
 *  compute RF = KF - (AF*UF + BF*U0)  *
 ***************************************/
mult_M24(&tempk1, &M->BF, &U->M[0]);
mult_tri2xu(&tempk2, &M->AF, &U->F);
addkv2(&tempk1, &tempk1, &tempk2);
subkv2(&R->F, &K->F, &tempk1);

/***********************************************
 *  compute R0 = K0 - (CF*UF + A0*U0 + B0*U1)  *
 ***********************************************/
mult_M42(&tempk1, &M->CF, &U->F);
mult_M24(&tempk2, &M->B[0], &U->M[1]);
seam(&k41, &tempk1, &tempk2);
mult_M44(&k42, &M->A[0], &U->M[0]);
addkv4(&k41, &k41, &k42);
subkv4(&R->M[0], &K->M[0], &k41);

/*******************************************************************
 *  compute R(i) = K(i) - (C(i-1)*U(i-1) + A(i)*U(i) + B(i)*U(i))  *
 *******************************************************************/
for (i = 1; i <= ym-3; i++)
    {
    im1 = i-1;
    mult_M24(&tempk1, &M->C[im1], &U->M[im1]);
    mult_M24(&tempk2, &M->B[i], &U->M[i+1]);
    seam(&k41, &tempk1, &tempk2);
    mult_M44(&k42, &M->A[i], &U->M[i]);
    addkv4(&k41, &k41, &k42);
    subkv4(&R->M[i], &K->M[i], &k41);
    }

/*********************************************************************************
 *  compute R(ym-2) = K(ym-2) - (C(ym-3)*U(ym-3) + A(ym-2)*U(ym-2) + B(L)*U(L))  *
 *********************************************************************************/
i = ym-2;
im1 = i-1;
mult_M24(&tempk1, &M->C[im1], &U->M[im1]);
mult_M42(&tempk2, &M->BL, &U->L);
seam(&k41, &tempk1, &tempk2);
mult_M44(&k42, &M->A[i], &U->M[i]);
addkv4(&k41, &k41, &k42);
subkv4(&R->M[i], &K->M[i], &k41);

/******************************************************
 *  compute R(L) = K(L) - (C(L)*U(ym-2) + A(L)*U(L))  *
 ******************************************************/
mult_M24(&tempk1, &M->CL, &U->M[ym-2]);
mult_tri2xu(&tempk2, &M->AL, &U->L);
addkv2(&tempk1, &tempk1, &tempk2);
subkv2(&R->L, &K->L, &tempk1);
}  /* end of function calc_res */


/**************************************************************
 *  FUNCTION:  norm2_collvec -- find the 2-norm of a collvec  *
 *  INPUT:  A = collvec                                       *
 *  OUTPUT: nm = norm                                         *
 **************************************************************/
double norm2_collvec (A)
struct collvec *A;
{
int i, j, k, l;
double nm = 0;

for (i = 0; i < xm; i++)
    for (j = 0; j < 2; j++)
        {
	nm += (A->F.el[i].el[j]) * (A->F.el[i].el[j]);
	nm += (A->L.el[i].el[j]) * (A->L.el[i].el[j]);
        }

for (i = 0; i < ym-1; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        nm += (A->M[i].el[j].el[k].el[l]) * (A->M[i].el[j].el[k].el[l]);

return(sqrt(nm));
}  /* end of function norm2_collvec */


/**************************************************************
 *  FUNCTION:  givsc -- change [x y]' by the Givens rotation  *
 *             defined by c and s                             *
 **************************************************************/
void givsc (x, y, c, s)
double *x, *y, c, s;
{
double tx = *x, ty = *y;

*x = c*tx - s*ty;
*y = s*tx + c*ty;
}  /* end of function givsc */


/*****************************************************************
 *  FUNCTION:  BCind  -- determine boundary unknowns and knowns  *
 *  INPUTS: unkNE, unkNW, unkSE, unkSW                           *
 *****************************************************************/
void BCind (unkNE, unkNW, unkSE, unkSW)
int *unkSW, *unkSE, *unkNW, *unkNE;
{
if (BCTypeW == 1)
    {
    if (BCTypeN == 1)
        *unkNW = 3;
    else
        *unkNW = 1;

    if (BCTypeS == 1)
        *unkSW = 3;
    else
        *unkSW = 1;
    }
else
    {
    if (BCTypeN == 1)
        *unkNW = 2;
    else
        *unkNW = 0;
    
    if (BCTypeS == 1)
        *unkSW = 2;
    else
        *unkSW = 0;
    }

if (BCTypeE == 1)
    {
    if (BCTypeN == 1)
        *unkNE = 3;
    else
        *unkNE = 1;

    if (BCTypeS == 1)
        *unkSE = 3;
    else
        *unkSE = 1;
    }
else
    {
    if (BCTypeN == 1)
        *unkNE = 2;
    else
        *unkNE = 0;

    if (BCTypeS == 1)
        *unkSE = 2;
    else
        *unkSE = 0;
    }
}  /* end of function BCind */


/***********************
 *  FUNCTION:  make_T  *
 ***********************/
void make_T (T, dx, dy, unkSW, unkSE, unkNW, unkNE, ccx, ccy, cx, cy)
struct collmat *T;
double *dx, *dy, **ccx, **ccy, **cx, **cy;
int unkSW, unkSE, unkNW, unkNE;
{
int i;

load_AF_CF_J(&T->AF, 0, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy);
load_AF_CF_J(&T->CF, 1, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy);

load_AL_BL_J(&T->AL, yc-1, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy);
load_AL_BL_J(&T->BL, yc-2, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy);

load_B_J(&T->BF, 0, 0, ccx, ccy, cx, cy, dx, dy);
for (i = 1; i <= ym-2; i++)
    load_B_J(&T->B[i-1], 2*i, i, ccx, ccy, cx, cy, dx, dy);

load_C_J(&T->CL, yc-1, ym-1, ccx, ccy, cx, cy, dx, dy);
for (i = 1; i <= ym-2; i++)
    load_C_J(&T->C[i-1], 2*i+1, i, ccx, ccy, cx, cy, dx, dy);

for (i = 1; i <= ym-1; i++)
    load_A_J(&T->A[i-1], 2*i, i, ccx, ccy, cx, cy, dx, dy);
}  /* end of function make_J */


/***********************
 *  FUNCTION:  make_M  *
 ***********************/
void make_M (M, dx, dy, unkSW, unkSE, unkNW, unkNE, ccx, ccy, cx, cy, t)
struct collmat *M;
double *dx, *dy, **ccx, **ccy, **cx, **cy, t;
int unkSW, unkSE, unkNW, unkNE;
{
int i;

load_AF_CF_L(&M->AF, 0, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy, t);
load_AF_CF_L(&M->CF, 1, unkSW, unkSE, ccx, ccy, cx, cy, dx, dy, t);

load_AL_BL_L(&M->AL, yc-1, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy, t);
load_AL_BL_L(&M->BL, yc-2, unkNW, unkNE, ccx, ccy, cx, cy, dx, dy, t);

load_B_L(&M->BF, 0, 0, ccx, ccy, cx, cy, dx, dy, t);
for (i = 1; i <= ym-2; i++)
    load_B_L(&M->B[i-1], 2*i, i, ccx, ccy, cx, cy, dx, dy, t);

load_C_L(&M->CL, yc-1, ym-1, ccx, ccy, cx, cy, dx, dy, t);
for (i = 1; i <= ym-2; i++)
    load_C_L(&M->C[i-1], 2*i+1, i, ccx, ccy, cx, cy, dx, dy, t);

for (i = 1; i <= ym-1; i++)
    load_A_L(&M->A[i-1], 2*i, i, ccx, ccy, cx, cy, dx, dy, t);
}  /* end of function make_M */


/***********************************************************************************************
 *  FUNCTION: load_BC_J -- load BC for diffrential operator J                                  *
 *  INPUTS:  SWc, SEc, NWc, NEc, Ws, Es, Ns, Ss = structures into which BC information goes    *
 *           dx, dy = vectors of element lengths                                               *
 ***********************************************************************************************/
void load_BC_J (SWc, SEc, NWc, NEc, Ws, Es, Ns, Ss, dx, dy, cx, cy, ccx, ccy)
struct corner_element_bcs *SWc, *SEc, *NWc, *NEc;
struct EW_element_bcs *Ws, *Es;
struct NS_element_bcs *Ns, *Ss;
double *dx, *dy, **cx, **cy, **ccx, **ccy;
{
int i, j, jp1, xind, yind, xsp, ysp;

/*  SW element  */
for (i = 0; i < 4; i++)
    {
    xind = i/2;
    yind = i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeS == 1)
        {
	SWc->el[i][0][0][0] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][0][0] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][0][0][1] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][0][1] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	SWc->el[i][0][0][0] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][0][0] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][0][0][1] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][0][1] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }

    if (BCTypeW == 1)
        {
	SWc->el[i][0][1][0] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][1][0] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][0][1][1] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][1][1] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	SWc->el[i][0][1][0] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][1][0] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][0][1][1] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SWc->el[i][1][1][1] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    }

/*  SE element  */
for (i = 0; i < 4; i++)
    {
    xind = xc - 2 + i/2;
    yind = i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeS == 1)
        {
	SEc->el[i][0][0][0] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][0][0] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][0][0][1] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][0][1] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	SEc->el[i][0][0][0] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][0][0] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][0][0][1] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][0][1] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }

    if (BCTypeE == 1)
        {
	SEc->el[i][0][1][0] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][1][0] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][0][1][1] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][1][1] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	SEc->el[i][0][1][0] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][1][0] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][0][1][1] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	SEc->el[i][1][1][1] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    }

/*  NW element  */
for (i = 0; i < 4; i++)
    {
    xind = i/2;
    yind = yc - 2 + i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeN == 1)
        {
	NWc->el[i][0][0][0] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][0][0] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][0][0][1] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][0][1] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	NWc->el[i][0][0][0] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][0][0] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][0][0][1] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][0][1] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    
    if (BCTypeW == 1)
        {
	NWc->el[i][0][1][0] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][1][0] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][0][1][1] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][1][1] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	NWc->el[i][0][1][0] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][1][0] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][0][1][1] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NWc->el[i][1][1][1] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    }

/*  NE element  */
for (i = 0; i < 4; i++)
    {
    xind = xc - 2 + i/2;
    yind = yc - 2 + i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeN == 1)
        {
	NEc->el[i][0][0][0] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][0][0] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][0][0][1] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][0][1] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	NEc->el[i][0][0][0] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][0][0] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][0][0][1] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][0][1] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }

    if (BCTypeE == 1)
        {
	NEc->el[i][0][1][0] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][1][0] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][0][1][1] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][1][1] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    else
        {
	NEc->el[i][0][1][0] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][1][0] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][0][1][1] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	NEc->el[i][1][1][1] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
        }
    }

/* W elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = i/2;
	yind = 2*jp1 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeW == 1)
	    {
	    Ws->el[j][i][0][0] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][1][0] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][0][1] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][1][1] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
	else
	    {
	    Ws->el[j][i][0][0] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][1][0] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][0][1] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ws->el[j][i][1][1] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
        }
    }

/* E elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = xc - 2 + i/2;
	yind = 2*jp1 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeE == 1)
	    {
	    Es->el[j][i][0][0] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][1][0] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][0][1] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][1][1] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
	else
	    {
	    Es->el[j][i][0][0] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][1][0] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][0][1] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Es->el[j][i][1][1] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
        }
    }

/* S elements */
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = 2*jp1 + i/2;
	yind = i%2;
	xsp = xind / 2;
	ysp = yind / 2;
	
	if (BCTypeS == 1)
	    {
	    Ss->el[j][i][0][0] = J_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][1][0] = J_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][0][1] = J_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][1][1] = J_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
	else
	    {
	    Ss->el[j][i][0][0] = J_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][1][0] = J_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][0][1] = J_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ss->el[j][i][1][1] = J_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
        }
    }

/* N elements */
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = 2*jp1 + i/2;
	yind = yc - 2 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeN == 1)
	    {
	    Ns->el[j][i][0][0] = J_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][1][0] = J_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][0][1] = J_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][1][1] = J_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
	else
	    {
	    Ns->el[j][i][0][0] = J_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][1][0] = J_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][0][1] = J_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    Ns->el[j][i][1][1] = J_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp]);
	    }
        }
    }
}  /* end of function load_BC_J */


/***********************************************************************************************
 *  FUNCTION: load_BC_L -- load BC for diffrential operator L                                  *
 *  INPUTS:  SWc, SEc, NWc, NEc, Ws, Es, Ns, Ss = structures into which BC information goes    *
 *           dx, dy = vectors of element lengths                                               *
 ***********************************************************************************************/
void load_BC_L (SWc, SEc, NWc, NEc, Ws, Es, Ns, Ss, dx, dy, cx, cy, ccx, ccy, t)
struct corner_element_bcs *SWc, *SEc, *NWc, *NEc;
struct EW_element_bcs *Ws, *Es;
struct NS_element_bcs *Ns, *Ss;
double *dx, *dy, **cx, **cy, **ccx, **ccy, t;
{
int i, j, jp1, xind, yind, xsp, ysp;

/*  SW element  */
for (i = 0; i < 4; i++)
    {
    xind = i/2;
    yind = i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeS == 1)
        {
	SWc->el[i][0][0][0] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][0][0] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][0][0][1] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][0][1] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	SWc->el[i][0][0][0] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][0][0] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][0][0][1] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][0][1] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }

    if (BCTypeW == 1)
        {
	SWc->el[i][0][1][0] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][1][0] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][0][1][1] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][1][1] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	SWc->el[i][0][1][0] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][1][0] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][0][1][1] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SWc->el[i][1][1][1] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    }

/*  SE element  */
for (i = 0; i < 4; i++)
    {
    xind = xc - 2 + i/2;
    yind = i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeS == 1)
        {
	SEc->el[i][0][0][0] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][0][0] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][0][0][1] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][0][1] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	SEc->el[i][0][0][0] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][0][0] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][0][0][1] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][0][1] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }

    if (BCTypeE == 1)
        {
	SEc->el[i][0][1][0] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][1][0] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][0][1][1] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][1][1] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	SEc->el[i][0][1][0] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][1][0] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][0][1][1] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	SEc->el[i][1][1][1] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    }

/*  NW element  */
for (i = 0; i < 4; i++)
    {
    xind = i/2;
    yind = yc - 2 + i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeN == 1)
        {
	NWc->el[i][0][0][0] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][0][0] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][0][0][1] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][0][1] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	NWc->el[i][0][0][0] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][0][0] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][0][0][1] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][0][1] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    
    if (BCTypeW == 1)
        {
	NWc->el[i][0][1][0] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][1][0] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][0][1][1] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][1][1] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	NWc->el[i][0][1][0] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][1][0] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][0][1][1] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NWc->el[i][1][1][1] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    }

/*  NE element  */
for (i = 0; i < 4; i++)
    {
    xind = xc - 2 + i/2;
    yind = yc - 2 + i%2;
    xsp = xind / 2;
    ysp = yind / 2;

    if (BCTypeN == 1)
        {
	NEc->el[i][0][0][0] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][0][0] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][0][0][1] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][0][1] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	NEc->el[i][0][0][0] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][0][0] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][0][0][1] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][0][1] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }

    if (BCTypeE == 1)
        {
	NEc->el[i][0][1][0] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][1][0] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][0][1][1] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][1][1] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    else
        {
	NEc->el[i][0][1][0] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][1][0] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][0][1][1] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	NEc->el[i][1][1][1] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
        }
    }

/* W elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = i/2;
	yind = 2*jp1 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeW == 1)
	    {
	    Ws->el[j][i][0][0] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][1][0] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][0][1] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][1][1] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
	else
	    {
	    Ws->el[j][i][0][0] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][1][0] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][0][1] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ws->el[j][i][1][1] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
        }
    }

/* E elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = xc - 2 + i/2;
	yind = 2*jp1 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeE == 1)
	    {
	    Es->el[j][i][0][0] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][1][0] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][0][1] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][1][1] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
	else
	    {
	    Es->el[j][i][0][0] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][1][0] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][0][1] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Es->el[j][i][1][1] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
        }
    }

/* S elements */
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = 2*jp1 + i/2;
	yind = i%2;
	xsp = xind / 2;
	ysp = yind / 2;
	
	if (BCTypeS == 1)
	    {
	    Ss->el[j][i][0][0] = L_fRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][1][0] = L_gRx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][0][1] = L_fLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][1][1] = L_gLx_fRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
	else
	    {
	    Ss->el[j][i][0][0] = L_fRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][1][0] = L_gRx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][0][1] = L_fLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ss->el[j][i][1][1] = L_gLx_gRy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
        }
    }

/* N elements */
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    for (i = 0; i < 4; i++)
        {
	xind = 2*jp1 + i/2;
	yind = yc - 2 + i%2;
	xsp = xind / 2;
	ysp = yind / 2;

	if (BCTypeN == 1)
	    {
	    Ns->el[j][i][0][0] = L_fRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][1][0] = L_gRx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][0][1] = L_fLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][1][1] = L_gLx_fLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
	else
	    {
	    Ns->el[j][i][0][0] = L_fRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][1][0] = L_gRx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][0][1] = L_fLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    Ns->el[j][i][1][1] = L_gLx_gLy(ccx[xind][yind], ccy[xind][yind], cx[xind][yind], cy[xind][yind], dx[xsp], dy[ysp] ,t);
	    }
        }
    }
}  /* end of function load_BC_L */

/*******************************************************************
 *  FUNCTION:  bdy_vals -- calculate boundary values               *
 *  INPUTS:  CWv, dCWv, CEv, dCEv, CSv, dCSv, CNv, dCNv = vectors  *
 *                       containing boundary values                *
 *           x, y = vectors containing node locations              *
 *******************************************************************/
void bdy_vals (Wv, dWv, Ev, dEv, Sv, dSv, Nv, dNv, x, y, t)
double *Wv, *dWv, *Ev, *dEv, *Sv, *dSv, *Nv, *dNv, *x, *y, t;
{
int i, j;

for (j = 0; j <= ym; j++)
    {
    Wv[j] =  BCValuW(y[j],t);
    dWv[j] = BCValudW(y[j],t);
    Ev[j] =  BCValuE(y[j],t);
    dEv[j] = BCValudE(y[j],t);
    }

for (i = 0; i <= xm; i++)
    {
    Sv[i] =  BCValuS(x[i],t);
    dSv[i] = BCValudS(x[i],t);
    Nv[i] =  BCValuN(x[i],t);
    dNv[i] = BCValudN(x[i],t);
    }
}  /*  end of function bdy_vals */


/****************************************************************
 *  FUNCTION:  sub_BC -- subtract BC's from forcing function    *
 *  INPUTS: HH = forcing function before and after subtraction  *
 ****************************************************************/
void sub_BC (HH, SWT, SET, NWT, NET, WT, ET, NT, ST, SWM, SEM, NWM, NEM, WM, EM, NM, SM, CWv, CEv, CNv, CSv, dCWv, dCEv, dCNv, dCSv, unkSW, unkSE, unkNW, unkNE, p, fact)
struct collvec *HH;
struct corner_element_bcs *SWM, *SEM, *NWM, *NEM, *SWT, *SET, *NWT, *NET;
struct EW_element_bcs *WM, *EM, *WT, *ET;
struct NS_element_bcs *NM, *SM, *NT, *ST;
double *CWv, *CEv, *CNv, *CSv, *dCWv, *dCEv, *dCNv, *dCSv, p, fact;
int unkSW, unkSE, unkNW, unkNE;
{
int i, j, jp1, jp2, yind;
double sigma[4];

/*  SW element  */
for (i = 0; i < 4; i++)
    {
    sigma[i] =   CWv[1] * (-p*SWM->el[i][0][1][0] + fact*SWT->el[i][0][1][0]);
    sigma[i] += dCWv[1] * (-p*SWM->el[i][1][1][0] + fact*SWT->el[i][1][1][0]);
    sigma[i] +=  CSv[1] * (-p*SWM->el[i][0][0][0] + fact*SWT->el[i][0][0][0]);
    sigma[i] += dCSv[1] * (-p*SWM->el[i][1][0][0] + fact*SWT->el[i][1][0][0]);
    sigma[i] +=  CSv[0] * (-p*SWM->el[i][0][0][1] + fact*SWT->el[i][0][0][1]);
    sigma[i] += dCSv[0] * (-p*SWM->el[i][1][0][1] + fact*SWT->el[i][1][0][1]);

    if (unkSW - BCTypeW == 2)
        sigma[i] += dCWv[0] * (-p*SWM->el[i][1][1][1] + fact*SWT->el[i][1][1][1]);
    else
        sigma[i] +=  CWv[0] * (-p*SWM->el[i][0][1][1] + fact*SWT->el[i][0][1][1]);
    }

HH->F.el[0].el[0]          += sigma[0];
HH->F.el[0].el[1]          += sigma[2];
HH->M[0].el[0].el[0].el[0] += sigma[1];
HH->M[0].el[0].el[1].el[0] += sigma[3];

/*  SE element  */
for (i = 0; i < 4; i++)
    {
    sigma[i] =   CEv[1]    * (-p*SEM->el[i][0][1][0] + fact*SET->el[i][0][1][0]);
    sigma[i] += dCEv[1]    * (-p*SEM->el[i][1][1][0] + fact*SET->el[i][1][1][0]);
    sigma[i] +=  CSv[xm-1] * (-p*SEM->el[i][0][0][0] + fact*SET->el[i][0][0][0]);
    sigma[i] += dCSv[xm-1] * (-p*SEM->el[i][1][0][0] + fact*SET->el[i][1][0][0]);
    sigma[i] +=  CSv[xm]   * (-p*SEM->el[i][0][0][1] + fact*SET->el[i][0][0][1]);
    sigma[i] += dCSv[xm]   * (-p*SEM->el[i][1][0][1] + fact*SET->el[i][1][0][1]);
    
    if (unkSE - BCTypeE == 2)
        sigma[i] += dCEv[0] * (-p*SEM->el[i][1][1][1] + fact*SET->el[i][1][1][1]);
    else
        sigma[i] +=  CEv[0] * (-p*SEM->el[i][0][1][1] + fact*SET->el[i][0][1][1]);
    }

HH->F   .el[xm-1].el[0]       += sigma[0];
HH->F   .el[xm-1].el[1]       += sigma[2];
HH->M[0].el[xm-1].el[0].el[0] += sigma[1];
HH->M[0].el[xm-1].el[1].el[0] += sigma[3];

/*  NW element  */
for (i = 0; i < 4; i++)
    {
    sigma[i] =   CWv[ym-1] * (-p*NWM->el[i][0][1][0] + fact*NWT->el[i][0][1][0]);
    sigma[i] += dCWv[ym-1] * (-p*NWM->el[i][1][1][0] + fact*NWT->el[i][1][1][0]);
    sigma[i] +=  CNv[1]    * (-p*NWM->el[i][0][0][0] + fact*NWT->el[i][0][0][0]);
    sigma[i] += dCNv[1]    * (-p*NWM->el[i][1][0][0] + fact*NWT->el[i][1][0][0]);
    sigma[i] +=  CNv[0]    * (-p*NWM->el[i][0][0][1] + fact*NWT->el[i][0][0][1]);
    sigma[i] += dCNv[0]    * (-p*NWM->el[i][1][0][1] + fact*NWT->el[i][1][0][1]);

    if (unkNW - BCTypeW == 2)
        sigma[i] += dCWv[ym] * (-p*NWM->el[i][1][1][1] + fact*NWT->el[i][1][1][1]);
    else
        sigma[i] +=  CWv[ym] * (-p*NWM->el[i][0][1][1] + fact*NWT->el[i][0][1][1]);
    }

yind = ym-2;
HH->L.el[0].el[0]             += sigma[1];
HH->L.el[0].el[1]             += sigma[3];
HH->M[yind].el[0].el[0].el[1] += sigma[0];
HH->M[yind].el[0].el[1].el[1] += sigma[2];

/*  NE element  */
for (i = 0; i < 4; i++)
    {
    sigma[i] =   CEv[ym-1] * (-p*NEM->el[i][0][1][1] + fact*NET->el[i][0][1][1]);
    sigma[i] += dCEv[ym-1] * (-p*NEM->el[i][1][1][1] + fact*NET->el[i][1][1][1]);
    sigma[i] +=  CNv[xm]   * (-p*NEM->el[i][0][0][0] + fact*NET->el[i][0][0][0]);
    sigma[i] += dCNv[xm]   * (-p*NEM->el[i][1][0][0] + fact*NET->el[i][1][0][0]);
    sigma[i] +=  CNv[xm-1] * (-p*NEM->el[i][0][0][1] + fact*NET->el[i][0][0][1]);
    sigma[i] += dCNv[xm-1] * (-p*NEM->el[i][1][0][1] + fact*NET->el[i][1][0][1]);
    
    if (unkNE - BCTypeE == 2)
        sigma[i] += dCEv[ym] * (-p*NEM->el[i][1][1][0] + fact*NET->el[i][1][1][0]);
    else
        sigma[i] +=  CEv[ym] * (-p*NEM->el[i][0][1][0] + fact*NET->el[i][0][1][0]);
    }

yind = ym-2;
HH->L      .el[xm-1].el[0]       += sigma[1];
HH->L      .el[xm-1].el[1]       += sigma[3];
HH->M[yind].el[xm-1].el[0].el[1] += sigma[0];
HH->M[yind].el[xm-1].el[1].el[1] += sigma[2];

/* W elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    jp2 = j+2;
    for (i = 0; i < 4; i++)
        {
	sigma[i] =   CWv[jp1] * (-p*WM->el[j][i][0][0] + fact*WT->el[j][i][0][0]);
	sigma[i] += dCWv[jp1] * (-p*WM->el[j][i][1][0] + fact*WT->el[j][i][1][0]);
	sigma[i] +=  CWv[jp2] * (-p*WM->el[j][i][0][1] + fact*WT->el[j][i][0][1]);
	sigma[i] += dCWv[jp2] * (-p*WM->el[j][i][1][1] + fact*WT->el[j][i][1][1]);
        }

    HH->M[j].el[0].el[0].el[1]   += sigma[0];
    HH->M[j].el[0].el[1].el[1]   += sigma[2];
    HH->M[jp1].el[0].el[0].el[0] += sigma[1];
    HH->M[jp1].el[0].el[1].el[0] += sigma[3];
    }

/* E elements */
for (j = 0; j <= ym-3; j++)
    {
    jp1 = j+1;
    jp2 = j+2;
    for (i = 0; i < 4; i++)
        {
	sigma[i] =   CEv[jp1] * (-p*EM->el[j][i][0][0] + fact*ET->el[j][i][0][0]);
	sigma[i] += dCEv[jp1] * (-p*EM->el[j][i][1][0] + fact*ET->el[j][i][1][0]);
	sigma[i] +=  CEv[jp2] * (-p*EM->el[j][i][0][1] + fact*ET->el[j][i][0][1]);
	sigma[i] += dCEv[jp2] * (-p*EM->el[j][i][1][1] + fact*ET->el[j][i][1][1]);
        }

    HH->M[j]  .el[xm-1].el[0].el[1] += sigma[0];
    HH->M[j]  .el[xm-1].el[1].el[1] += sigma[2];
    HH->M[jp1].el[xm-1].el[0].el[0] += sigma[1];
    HH->M[jp1].el[xm-1].el[1].el[0] += sigma[3];
    }

/* S elements */
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    jp2 = j+2;
    for (i = 0; i < 4; i++)
        {
	sigma[i] =   CSv[jp1] * (-p*SM->el[j][i][0][0] + fact*ST->el[j][i][0][0]);
	sigma[i] += dCSv[jp1] * (-p*SM->el[j][i][1][0] + fact*ST->el[j][i][1][0]);
	sigma[i] +=  CSv[jp2] * (-p*SM->el[j][i][0][1] + fact*ST->el[j][i][0][1]);
	sigma[i] += dCSv[jp2] * (-p*SM->el[j][i][1][1] + fact*ST->el[j][i][1][1]);
        }

    HH->F.   el[jp1].el[0]       += sigma[0];
    HH->F.   el[jp1].el[1]       += sigma[2];
    HH->M[0].el[jp1].el[0].el[0] += sigma[1];
    HH->M[0].el[jp1].el[1].el[0] += sigma[3];
    }
  
/* N elements */
yind = ym-2;
for (j = 0; j <= xm-3; j++)
    {
    jp1 = j+1;
    jp2 = j+2;
    for (i = 0; i < 4; i++)
        {
	sigma[i] =   CNv[jp1] * (-p*NM->el[j][i][0][0] + fact*NT->el[j][i][0][0]);
	sigma[i] += dCNv[jp1] * (-p*NM->el[j][i][1][0] + fact*NT->el[j][i][1][0]);
	sigma[i] +=  CNv[jp2] * (-p*NM->el[j][i][0][1] + fact*NT->el[j][i][0][1]);
	sigma[i] += dCNv[jp2] * (-p*NM->el[j][i][1][1] + fact*NT->el[j][i][1][1]);
        }

    HH->L      .el[jp1].el[0]       += sigma[1];
    HH->L      .el[jp1].el[1]       += sigma[3];
    HH->M[yind].el[jp1].el[0].el[1] += sigma[0];
    HH->M[yind].el[jp1].el[1].el[1] += sigma[2];
    }
}  /* end of function sub_BC */


/*******************************************************************
 *  FUNCTION:  zap_collvec -- sets every element of a collvec = 0  *
 *  INPUT: A = collvec                                             *
 *******************************************************************/
void zap_collvec (A)
struct collvec *A;
{
int i, j, k, l;

for (i = 0; i < xm; i++)
    for (j = 0; j < 2; j++)
        {
	A->F.el[i].el[j] = 0;
	A->L.el[i].el[j] = 0;
        }

for (i = 0; i <= ym-2; i++)
    for (j = 0; j < xm; j++)
        for (k = 0; k < 2; k++)
	    for (l = 0; l < 2; l++)
	        A->M[i].el[j].el[k].el[l] = 0;
}  /* end of function zap_collvec  */


/**************************************************************
 *  FUNCTION: pbcgstab -- solves Ax=b by preconditioned       *
 *                   bi-conjugate gradient method stabilized  *
 **************************************************************/
void pbcgstab (x, A, b, iter, nm)
struct collvec *x, *b;
struct collmat *A;
int *iter;
double *nm;
{
struct collvec r, rbar, v, p, y, s, t, z;
double rho, oldrho, alpha, omega, beta;

mult_collvec_collmat(&r, A, x);
sub_collvec(&r, b, &r);
copy_collvec(&rbar, &r);

rho = 1;
alpha = 1;
omega = 1;

zap_collvec(&v);
zap_collvec(&p);

*iter = 0;
while (*iter <= maxiter)
    {
    (*iter)++;

    oldrho = rho;
    rho = dot_collvec(&rbar, &r);
    beta = (rho/oldrho) * (alpha/omega);
    
    saxpy_collvec(&p, -omega, &v, &p);
    saxpy_collvec(&p, beta, &p, &r);
    
    psolve(&y, A, &p);
    mult_collvec_collmat_trick(&v, A, &y, &p);

    alpha = rho / dot_collvec(&rbar, &v);
    saxpy_collvec(&s, -alpha, &v, &r);
    
    psolve(&z, A, &s);    
    mult_collvec_collmat(&t, A, &z, &s);

    omega = dot_collvec(&t, &s) / dot_collvec(&t, &t);
    
    saxpy_collvec(x, alpha, &y, x);
    saxpy_collvec(x, omega, &z, x);
    
    saxpy_collvec(&r, -omega, &t, &s);
    
    *nm = norm_collvec(&r);
    if (pr)
        printf("%d   %e\n", *iter, *nm);
    
    if (*nm <= restol)
        break;
    }
}  /* end of function pbcgstab */


/************************************************
 *  FUNCTION: gmres -- solves Ax=b by GMRES(k)  *
 ************************************************/
void gmres (x, A, b, iter, nm)
struct collvec *x, *b;
struct collmat *A;
int *iter;
double *nm;
{
int conv = FALSE, i, im1, j, p, ctr, ix, reps;
double *Hvec, *be1, beta, c, s, t, *cc, *ss, aa, bb, *y;
struct collvec tmpcv, V[restart+1], r, w, sm, vhat;

/**************
 *  allocate  *
 **************/
cc = dvector(restart);
ss = dvector(restart);
y = dvector(restart);

reps = 0;
while (!conv && *iter <= maxiter)
    {
    reps++;
    blank(1);

    /*********
     *  zap  *
     *********/
    for (i = 0; i <= restart; i++)
        zap_collvec(&V[i]);
    Hvec = dvector(restart*(restart+3)/2);
    be1 = dvector(restart+1);

    /****************
     *  initialize  *
     ****************/
    calc_res(&tmpcv, A, x, b);
    psolve(&r, A, &tmpcv);
    beta = norm2_collvec(&r);
    scmult_collvec (&V[0], &r, 1/beta);
    be1[0] = beta;

    j = 0;
    ctr = -1;
    while ((j < restart) && (!conv))
        {
	(*iter)++;

	/****************************************************
	 *  Arnoldi to generate V with orthonormal columns  *
	 ****************************************************/
	mult_collvec_collmat(&tmpcv, A, &V[j]);
	psolve(&w, A, &tmpcv);

	zap_collvec(&sm);
	for (i = 0; i <= j; i++)
	    {
	    ctr++;
	    Hvec[ctr] = dot_collvec(&w, &V[i]);
	    saxpy_collvec(&sm, Hvec[ctr], &V[i], &sm);
	    
	    if (i > 0)
	        {
		im1 = i-1;
		givsc(&Hvec[ctr-1], &Hvec[ctr], cc[im1], ss[im1]);
	        }
	    }
	
	sub_collvec(&vhat, &w, &sm);
	ctr++;
	Hvec[ctr] = norm2_collvec(&vhat);
	scmult_collvec (&V[j+1], &vhat, 1/Hvec[ctr]);

	/**************************************************************
	 *  Givens rotation: H: upper Hessenberg -> upper triangular  *
	 **************************************************************/
	ix = j * (j+5) / 2;
	aa = Hvec[ix];
	bb = Hvec[ix+1];

	if (bb == 0)
	    {
	    c = 1;
	    s = 0;
	    }
	else
	    {
	    if (fabs(bb) > fabs(aa))
	        {
		t = -aa/bb;
		s = 1 / sqrt(1+t*t);
		c = s * t;
	        }
	    else
	        {
		t = -bb/aa;
		c = 1 / sqrt(1+t*t);
		s = c * t;
	        }
	    }
	
	Hvec[ix] = c*Hvec[ix] - s*Hvec[ix+1];
	Hvec[ix+1] = 0;

	givsc(&be1[j], &be1[j+1], c, s);
	
	cc[j] = c;
	ss[j] = s;
	
	*nm = fabs(be1[j+1]);
	if (pr) 
	    printf("reps = %d, j = %d, norm = %e\n", reps, j, *nm);
	
	if (*nm <= restol)
	    conv = TRUE;
	
	j++;
        }

    /***********************************************
     *  back substitution and compute new vector x *
     ***********************************************/
    for (i = j-1; i >= 0; i--)
        {
	y[i] = be1[i];
	for (p = i+1; p <= j; p++)
	    y[i] -= Hvec[i+p*(p+3)/2] * y[p];
	y[i] /= Hvec[i*(i+5)/2];
	saxpy_collvec(x, y[i], &V[i], x);
        }

    /*****************
     *  de-allocate  *
     *****************/
    free(Hvec);
    free(be1);
    }

free(cc);
free(ss);
free(y);
}  /* end of function gmres */


/***************************************************************************
 *  FUNCTION:  soln_gen -- converts the solution vector u and BC's         *
 *                              to a pleasant format:                      *
 *                                                                         *
 *                  time = ________                                        *
 *                  u =                                                    *
 *                  u_x =                                                  *
 *                  u_y =                                                  *
 *                  u_xy =                                                 *
 *  INPUTS:  U = solution vector for unknowns                              *
 *           CWv, dCWv, CEv, dCEv, CNv, dCNv,CSv, dCSv, = vectors of BC's  *
 *           S = 3-dimensional matrix containing solution                  *
 ***************************************************************************/
void soln_gen (S, U, CWv, dCWv, CEv, dCEv, CSv, dCSv, CNv, dCNv, unkSW, unkSE, unkNW, unkNE)
struct collvec *U;
double ***S, *CWv, *dCWv, *CEv, *dCEv, *CSv, *dCSv, *CNv, *dCNv;
int unkSW, unkSE, unkNW, unkNE;
{
int i, j, k;

/*  interior nodes  */
for (j = 1; j < ym; j++)
    {
    k = ym-j-1;
    for (i = 1; i < xm; i++)
        {
	S[0][j][i] = U->M[k].el[i-1].el[1].el[0];
	S[1][j][i] = U->M[k].el[i]  .el[0].el[0];
	S[2][j][i] = U->M[k].el[i-1].el[1].el[1];
	S[3][j][i] = U->M[k].el[i]  .el[0].el[1];
        }
    }

/*  E nodes  */
for (j = 1; j < ym; j++)
    {
    k = ym-j;
    S[1-BCTypeE][j][xm] =  CEv[k];
    S[3-BCTypeE][j][xm] = dCEv[k];
    }
for (j = 1; j < ym; j++)
    {
    k = ym-j-1;
    S[BCTypeE]  [j][xm] = U->M[k].el[xm-1].el[1].el[0];
    S[BCTypeE+2][j][xm] = U->M[k].el[xm-1].el[1].el[1];
    }

/*  W nodes  */
for (j = 1; j < ym; j++)
    {
    k = ym-j;
    S[1-BCTypeW]  [j][0] =  CWv[k];
    S[3-BCTypeW][j][0] = dCWv[k];
    }
for (j = 1; j < ym; j++)
    {
    k = ym-j-1;
    S[BCTypeW]  [j][0] = U->M[k].el[0].el[0].el[0];
    S[BCTypeW+2][j][0] = U->M[k].el[0].el[0].el[1];
    }

/*  S nodes  */
k = 2*(1-BCTypeS);
for (i = 1; i < xm; i++)
    {
    S[k]  [ym][i] =  CSv[i];
    S[k+1][ym][i] = dCSv[i];
    }
k = 2*BCTypeS;
for (i = 1; i < xm; i++)
    {
    S[k]  [ym][i] = U->F.el[i-1].el[1];
    S[k+1][ym][i] = U->F.el[i]  .el[0];
    }

/*  N nodes  */
k = 2*(1-BCTypeN);
for (i = 1; i < xm; i++)
    {
    S[k]  [0][i] =  CNv[i];
    S[k+1][0][i] = dCNv[i];
    }
k = 2*BCTypeN;
for (i = 1; i < xm; i++)
    {
    S[k]  [0][i] = U->L.el[i-1].el[1];
    S[k+1][0][i] = U->L.el[i]  .el[0];
    }

/*  SW corner node  */
if (unkSW == 0)
    {
    S[0][ym][0] = U->F.el[0].el[0];
    S[1][ym][0] =  CWv[0];
    S[2][ym][0] =  CSv[0];
    S[3][ym][0] = dCSv[0];
    }
else if (unkSW == 1)
    {
    S[0][ym][0] =  CWv[0];
    S[1][ym][0] = U->F.el[0].el[0];
    S[2][ym][0] =  CSv[0];
    S[3][ym][0] = dCSv[0];
    }
else if (unkSW == 2)
    {
    S[0][ym][0] =  CSv[0];
    S[1][ym][0] = dCSv[0];
    S[2][ym][0] = U->F.el[0].el[0];
    S[3][ym][0] = dCWv[0];
    }
else
    {
    S[0][ym][0] =  CSv[0];
    S[1][ym][0] = dCSv[0];
    S[2][ym][0] = dCWv[0];
    S[3][ym][0] = U->F.el[0].el[0];
    }

/*  SE corner node  */
if (unkSE == 0)
    {
    S[0][ym][xm] = U->F.el[xm-1].el[1];
    S[1][ym][xm] =  CEv[0];
    S[2][ym][xm] =  CSv[xm];
    S[3][ym][xm] = dCSv[xm];
    }
else if (unkSE == 1)
    {
    S[0][ym][xm] =  CEv[0];
    S[1][ym][xm] = U->F.el[xm-1].el[1];
    S[2][ym][xm] =  CSv[xm];
    S[3][ym][xm] = dCSv[xm];
    }
else if (unkSE == 2)
    {
    S[0][ym][xm] =  CSv[xm];
    S[1][ym][xm] = dCSv[xm];
    S[2][ym][xm] = U->F.el[xm-1].el[1];
    S[3][ym][xm] = dCEv[0];
    }
else
    {
    S[0][ym][xm] =  CSv[xm];
    S[1][ym][xm] = dCSv[xm];
    S[2][ym][xm] = dCEv[0];
    S[3][ym][xm] = U->F.el[xm-1].el[1];
    }

/*  NW corner node  */
if (unkNW == 0)
    {
    S[0][0][0] = U->L.el[0].el[0];
    S[1][0][0] =  CWv[ym];
    S[2][0][0] =  CNv[0];
    S[3][0][0] = dCNv[0];
    }
else if (unkNW == 1)
    {
    S[0][0][0] =  CWv[ym];
    S[1][0][0] = U->L.el[0].el[0];
    S[2][0][0] =  CNv[0];
    S[3][0][0] = dCNv[0];
    }
else if (unkNW == 2)
    {
    S[0][0][0] =  CNv[0];
    S[1][0][0] = dCNv[0];
    S[2][0][0] = U->L.el[0].el[0];
    S[3][0][0] = dCWv[ym];
    }
else
    {
    S[0][0][0] =  CNv[0];
    S[1][0][0] = dCNv[0];
    S[2][0][0] = dCWv[ym];
    S[3][0][0] = U->L.el[0].el[0];
    }

/*  NE corner node  */
if (unkNE == 0)
    {
    S[0][0][xm] = U->L.el[xm-1].el[1];
    S[1][0][xm] =  CEv[ym];
    S[2][0][xm] =  CNv[xm];
    S[3][0][xm] = dCNv[xm];
    }
else if (unkNE == 1)
    {
    S[0][0][xm] =  CEv[ym];
    S[1][0][xm] = U->L.el[xm-1].el[1];
    S[2][0][xm] =  CNv[xm];
    S[3][0][xm] = dCNv[xm];
    }
else if (unkNE == 2)
    {
    S[0][0][xm] =  CNv[xm];
    S[1][0][xm] = dCNv[xm];
    S[2][0][xm] = U->L.el[xm-1].el[1];
    S[3][0][xm] = dCEv[ym];
    }
else
    {
    S[0][0][xm] =  CNv[xm];
    S[1][0][xm] = dCNv[xm];
    S[2][0][xm] = dCEv[ym];
    S[3][0][xm] = U->L.el[xm-1].el[1];
    }
}  /* end of fucntion soln_gen */


/**************************************************************************
 *  FUNCTION:  load_collvec_ic -- load a collvec with initial conditions  *
 **************************************************************************/
void load_collvec_ic (U, x, y, unkSW, unkSE, unkNW, unkNE, t)
struct collvec *U;
double *x, *y, t;
int unkSW, unkSE, unkNW, unkNE;
{
int i, j;

for (i = 1; i <= ym-1; i++)
    {
    for (j = 1; j <= xm-1; j++)
        {
	U->M[i-1].el[j-1].el[1].el[0] =    S(x[j], y[i], t);
	U->M[i-1].el[j-1].el[1].el[1] =  S_y(x[j], y[i], t);
	U->M[i-1].el[j]  .el[0].el[0] =  S_x(x[j], y[i], t);
	U->M[i-1].el[j]  .el[0].el[1] = S_xy(x[j], y[i], t);
        }

    if (BCTypeW == 1)
        {
	U->M[i-1].el[0].el[0].el[0] =  S_x(x[0], y[i], t);
	U->M[i-1].el[0].el[0].el[1] = S_xy(x[0], y[i], t);
        }
    else
        {
	U->M[i-1].el[0].el[0].el[0] =   S(x[0], y[i], t);
	U->M[i-1].el[0].el[0].el[1] = S_y(x[0], y[i], t);
        }

    if (BCTypeE == 1)
        {
	U->M[i-1].el[xm-1].el[1].el[0] =  S_x(x[xm], y[i], t);
	U->M[i-1].el[xm-1].el[1].el[1] = S_xy(x[xm], y[i], t);
        }
    else
        {
	U->M[i-1].el[xm-1].el[1].el[0] =   S(x[xm], y[i], t);
	U->M[i-1].el[xm-1].el[1].el[1] = S_y(x[xm], y[i], t);
        }
    }

for (j = 1; j <= xm-1; j++)
    {
    if (BCTypeS == 1)
        {
	U->F.el[j-1].el[1] =  S_y(x[j], y[0], t);
	U->F.el[j]  .el[0] = S_xy(x[j], y[0], t);
        }
    else
        {
	U->F.el[j-1].el[1] =   S(x[j], y[0], t);
	U->F.el[j]  .el[0] = S_x(x[j], y[0], t);
        }

    if (BCTypeN == 1)
        {
	U->L.el[j-1].el[1] =  S_y(x[j], y[ym], t);
	U->L.el[j]  .el[0] = S_xy(x[j], y[ym], t);
        }
    else
        {
	U->L.el[j-1].el[1] =   S(x[j], y[ym], t);
	U->L.el[j]  .el[0] = S_x(x[j], y[ym], t);
        }
    }

if (unkSW == 0)
    U->F.el[0].el[0] =         S(x[0], y[0], t);
else if (unkSW == 1)
    U->F.el[0].el[0] =       S_x(x[0], y[0], t);
else if (unkSW == 2)
    U->F.el[0].el[0] =       S_y(x[0], y[0], t);
else
    U->F.el[0].el[0] =      S_xy(x[0], y[0], t);

if (unkSE == 0)
    U->F.el[xm-1].el[1] =     S(x[xm], y[0], t);
else if (unkSE == 1)
    U->F.el[xm-1].el[1] =   S_x(x[xm], y[0], t);
else if (unkSE == 2)
    U->F.el[xm-1].el[1] =   S_y(x[xm], y[0], t);
else
    U->F.el[xm-1].el[1] =  S_xy(x[xm], y[0], t);

if (unkNW == 0)
    U->L.el[0].el[0] =        S(x[0], y[ym], t);
else if (unkNW == 1)
    U->L.el[0].el[0] =      S_x(x[0], y[ym], t);
else if (unkNW == 2)
    U->L.el[0].el[0] =      S_y(x[0], y[ym], t);
else
    U->L.el[0].el[0] =     S_xy(x[0], y[ym], t);

if (unkNE == 0)
    U->L.el[xm-1].el[1] =    S(x[xm], y[ym], t);
else if (unkNE == 1)
    U->L.el[xm-1].el[1] =  S_x(x[xm], y[ym], t);
else if (unkNE == 2)
    U->L.el[xm-1].el[1] =  S_y(x[xm], y[ym], t);
else
    U->L.el[xm-1].el[1] = S_xy(x[xm], y[ym], t);


}  /* end of fucntion load_collvec_ic */


/********************************
 *  FUNCTION: saxpy_triblock22  *
 ********************************/
void saxpy_triblock22 (Z, a, X, Y)
struct triblock22 *Z, *X, *Y;
double a;
{
int i, j, k;

for (k = 0; k < xm; k++)
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    Z->a[k].el[i][j] = a*X->a[k].el[i][j] + Y->a[k].el[i][j];

for (k = 0; k < xm-1; k++)
    for (j = 0; j < 2; j++)
        {
	Z->b[k].el[j] = a*X->b[k].el[j] + Y->b[k].el[j];
	Z->c[k].el[j] = a*X->c[k].el[j] + Y->c[k].el[j];
        }
}  /*  end of function saxpy_triblock22 */


/********************************
 *  FUNCTION: saxpy_triblock24  *
 ********************************/
void saxpy_triblock24 (Z, a, X, Y)
struct triblock24 *Z, *X, *Y;
double a;
{
int i, j, k, p;

for (k = 0; k < xm; k++)
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    for (p = 0; p < 2; p++)
	        Z->a[k].el[p].el[i][j] = a*X->a[k].el[p].el[i][j] + Y->a[k].el[p].el[i][j];

for (k = 0; k < xm-1; k++)
    for (j = 0; j < 2; j++)
        for (p = 0; p < 2; p++)
	    {
	    Z->b[k].el[j][p] = a*X->b[k].el[j][p] + Y->b[k].el[j][p];
	    Z->c[k].el[j][p] = a*X->c[k].el[j][p] + Y->c[k].el[j][p];
	    }
}  /*  end of function saxpy_triblock24 */


/********************************
 *  FUNCTION: saxpy_triblock44  *
 ********************************/
void saxpy_triblock44 (Z, a, X, Y)
struct triblock44 *Z, *X, *Y;
double a;
{
int i, j, k, p;

for (k = 0; k < xm; k++)
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    {
	    Z->a[k].A.el[i][j] = a*X->a[k].A.el[i][j] + Y->a[k].A.el[i][j];
	    Z->a[k].B.el[i][j] = a*X->a[k].B.el[i][j] + Y->a[k].B.el[i][j];
	    Z->a[k].C.el[i][j] = a*X->a[k].C.el[i][j] + Y->a[k].C.el[i][j];
	    Z->a[k].D.el[i][j] = a*X->a[k].D.el[i][j] + Y->a[k].D.el[i][j];
	    }

for (k = 0; k < xm-1; k++)
    for (j = 0; j < 2; j++)
        for (i = 0; i < 2; i++)
	    for (p = 0; p < 2; p++)
	    {
	    Z->b[k].el[j].el[i][p] = a*X->b[k].el[j].el[i][p] + Y->b[k].el[j].el[i][p];
	    Z->c[k].el[j].el[i][p] = a*X->c[k].el[j].el[i][p] + Y->c[k].el[j].el[i][p];
	    }
}  /*  end of function saxpy_triblock44 */



/******************************
 *  FUNCTION:  saxpy_collmat  *
 ******************************/
saxpy_collmat(G, a, M, T)
struct collmat *G, *M, *T;
double a;
{
int i;

saxpy_triblock22(&G->AF, a, &M->AF, &T->AF);
saxpy_triblock22(&G->CF, a, &M->CF, &T->CF);
saxpy_triblock22(&G->AL, a, &M->AL, &T->AL);
saxpy_triblock22(&G->BL, a, &M->BL, &T->BL);

saxpy_triblock24(&G->BF, a, &M->BF, &T->BF);
saxpy_triblock24(&G->CL, a, &M->CL, &T->CL);
for (i = 0; i <= ym-3; i++)
    {
    saxpy_triblock24(&G->B[i], a, &M->B[i], &T->B[i]);
    saxpy_triblock24(&G->C[i], a, &M->C[i], &T->C[i]);
    }

for (i = 0; i <= ym-2; i++)
    saxpy_triblock44(&G->A[i], a, &M->A[i], &T->A[i]);
}  /*  end of function saxpy_collmat */