/*************************************************************************************
 *  FUNCTION:  file2split_collvec --  read entries for a collvec from a file, place  *
 *         them in vectors where each vector will be sent to a particular processor  *
 *************************************************************************************/
void file2split_collvec (V, filename, nproc, bpp, xm, ym, num, loc, start, size)
double *V;
char filename[80];
int nproc, bpp, xm, ym, num, loc, start, size;
{
int i, j, k, l, p, stp, *ctr;
FILE *fp;

ctr = ivector(nproc);
for (i = 0; i < nproc; i++)
    ctr[i] = start - 1;

fp = fopen(filename, "r");

if ((xm<ym && num==0) || (xm>=ym && num==1))  /* v2cvA */
    {
    /***   FIRST   ***/
    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &V[++(ctr[0])]);

    /***   MIDDLE   ***/
    for (p = 0; p < nproc; p++)
        {
	if (p == 0)
	    stp = bpp-2;
	else
	    stp = bpp-1;
        for (i = 0; i <= stp; i++)
	    for (l = 0; l < 2; l++)
	        for (j = 0; j < loc; j++)
		    for (k = 0; k < 2; k++)
		        fscanf(fp, "%lf", &V[size*p+(++(ctr[p]))]);
        }

    /***   LAST   ***/
    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &V[++(ctr[0])]);
    }
else                                          /* v2cvB */
    {
    for (l = 0; l < loc; l++)
        for (i = 0; i < 2; i++)
	    {
	    /***   FIRST   ***/
	    fscanf(fp, "%lf", &V[++(ctr[0])]);

	    /***   MIDDLE   ***/
	    for (p = 0; p < nproc; p++)
	        {
		if (p == 0)
		    stp = bpp-2;
		else
		    stp = bpp-1;
		for (j = 0; j <= stp; j++)
		    for (k = 0; k < 2; k++)
		        fscanf(fp, "%lf", &V[size*p+(++(ctr[p]))]);
	        }

	    /***   LAST   ***/
	    fscanf(fp, "%lf", &V[++(ctr[0])]);
	    }
    }

fclose(fp);

free(ctr);
}  /* end of function file2_splitcollvec */


/*******************************
 *  FUNCTION:  load_collvec_M  *
 *******************************/
void load_collvec_M (M, V, bpp, xm, ym, num, loc, start)
struct collvec *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int i, j, k, l, ctr;

ctr = start - 1;

if ((xm<ym && num==0) || (xm>=ym && num==1))  /* v2cvA */
    {
    /***   FIRST   ***/
    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    M->F[0].el[i].el[j] = V[++ctr];

    /***   MIDDLE   ***/
    for (i = 0; i <= bpp-2; i++)
        for (l = 0; l < 2; l++)
	    for (j = 0; j < loc; j++)
	        for (k = 0; k < 2; k++)
		    M->M[i].el[j].el[k].el[l] = V[++ctr];

    /***   LAST   ***/
    for (i = 0; i < loc; i++)
        for (j = 0; j < 2; j++)
	    M->L[0].el[i].el[j] = V[++ctr];
    }
else                                          /* v2cvB */
    {
    for (l = 0; l < loc; l++)
        for (i = 0; i < 2; i++)
	    {
	    /***   FIRST   ***/
	    M->F[0].el[l].el[i] = V[++ctr];

	    /***   MIDDLE   ***/
	    for (j = 0; j <= bpp-2; j++)
	        for (k = 0; k < 2; k++)
		    M->M[j].el[l].el[i].el[k] = V[++(ctr)];

	    /***   LAST   ***/
	    M->L[0].el[l].el[i] = V[++ctr];
	    }
    }

}  /* end of function load_collvec_M */


/*******************************
 *  FUNCTION:  load_collvec_S  *
 *******************************/
void load_collvec_S (M, V, bpp, xm, ym, num, loc, start)
struct mat_kvec4 *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int i, j, k, l, ctr;

ctr = start - 1;

if ((xm<ym && num==0) || (xm>=ym && num==1))  /* v2cvA */
    {
    for (i = 0; i <= bpp-1; i++)
        for (l = 0; l < 2; l++)
	    for (j = 0; j < loc; j++)
	        for (k = 0; k < 2; k++)
		    M->el[i].el[j].el[k].el[l] = V[++ctr];
    }
else                                          /* v2cvB */
    {
    for (l = 0; l < loc; l++)
        for (i = 0; i < 2; i++)
	    for (j = 0; j <= bpp-1; j++)
	        for (k = 0; k < 2; k++)
		    M->el[j].el[l].el[i].el[k] = V[++(ctr)];
    }

}  /* end of function load_collvec_S */


/*********************************
 *  FUNCTION:  unload_collvec_S  *
 *********************************/
void unload_collvec_S (V, M, bpp, xm, ym, num, loc, start)
struct mat_kvec4 *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int i, j, k, l, ctr;

ctr = start - 1;

if ((xm<ym && num==0) || (xm>=ym && num==1))  /* v2cvA */
    {
    for (i = 0; i <= bpp-1; i++)
        for (l = 0; l < 2; l++)
	    for (j = 0; j < loc; j++)
	        for (k = 0; k < 2; k++)
		    V[++ctr] = M->el[i].el[j].el[k].el[l];
    }
else                                          /* v2cvB */
    {
    for (l = 0; l < loc; l++)
        for (i = 0; i < 2; i++)
	    for (j = 0; j <= bpp-1; j++)
	        for (k = 0; k < 2; k++)
		    V[++ctr] = M->el[j].el[l].el[i].el[k];
    }

}  /* end of function unload_collvec_S */


/*************************************************************
 *  FUNCTION:  load_SW_el_split -- loads the SW element      *
 *              of a collmat into vectors where each vector  *
 *              will be sent to a particular processor       *
 *  M = collmat                                              *
 *  f = 0, use xl_yg numbering                               *
 *  f = 1, use yl_xg numbering                               *
 *  ctr = keeps track for each processor                     *
 *************************************************************/
void load_SW_el_split (M, fp, f, ctr, size)
double *M;
FILE *fp;
int f, *ctr, size;
{
int i, j, k;

if (f == 0)  /* xl_yg */
    {
    /*  SW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  NW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  NE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    }
else  /* f == 1, yl_xg */
    {
    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
}  /* end of function load_SW_el_split */


/********************************************************************
 *  FUNCTION:  load_S_el_split -- loads the S element of a collmat  *
 *  M = collmat                                                     *
 *  f = 0, use xl_yg numbering                                      *
 *  f = 1, use yl_xg numbering                                      *
 *  m = index for which S element                                   *
 ********************************************************************/
void load_S_el_split (M, fp, f, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, m, nproc, bpp, *ctr, size;
{
int i, j, ix, ixp1;

if (f == 0)  /* xl_yg */
    {
    /*  SW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  NW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  NE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    }
else  /* f == 1, yl_xg */
    {
    ix   = (m+1) / bpp;
    ixp1 = (m+2) / bpp;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);

    /*  SE equation  */
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);

    /*  NE equation  */
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    }
}  /* end of function load_S_el_split */


/**********************************************************************
 *  FUNCTION:  load_SE_el_split -- loads the SE element of a collmat  *
 *  M = collmat                                                       *
 *  f = 0, use xl_yg numbering                                        *
 *  f = 1, use yl_xg numbering                                        *
 *  m = index for which S element                                     *
 **********************************************************************/
void load_SE_el_split (M, fp, f, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, m, nproc, bpp, *ctr, size;
{
int i, j, k, ix, nslaves;

if (f == 0)  /* xl_yg */
    {
    /*  SW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    
    /*  NW equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  SE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
else  /* f == 1, yl_xg */
    {
    ix = (m+1) / bpp;
    nslaves = nproc-1;
      
    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    }
}  /* end of function load_SE_el_split */


/********************************************************************
 *  FUNCTION:  load_W_el_split -- loads the W element of a collmat  *
 *  M = collmat                                                     *
 *  f = 0, use xl_yg numbering                                      *
 *  f = 1, use yl_xg numbering                                      *
 *  n = index for which W element                                   *
 ********************************************************************/
void load_W_el_split (M, fp, f, n, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, nproc, bpp, *ctr, size;
{
int i, j, k, nm1, ix, ixp1;

if (f == 0)  /* xl_yg */
    {
    ix   = (n+1) / bpp;
    ixp1 = (n+2) / bpp;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
	fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
        }

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
	fscanf(fp, "%lf", &M[size*ixp1+(++(ctr[ixp1]))]);
        }
    }
else  /* f == 1, yl_xg */
    {
    nm1 = n-1;

    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
}  /* end of fucntion load_W_el_split */


/****************************************************************************
 *  FUNCTION:  load_int_el_split -- loads an interior element of a collmat  *
 *  M = collmat                                                             *
 *  f = 0, use xl_yg numbering                                              *
 *  f = 1, use yl_xg numbering                                              *
 *  n, m = indices for which element                                        *
 ****************************************************************************/
void load_int_el_split (M, fp, f, n, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, m, nproc, bpp, *ctr, size;
{
int i, j, nx, nxp1, mx, mxp1;

if (f == 0)  /* xl_yg */
    {
    nx   = (n+1) / bpp;
    nxp1 = (n+2) / bpp;
    
    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }
    }
else  /* f == 1, yl_xg */
    {
    mx   = (m+1) / bpp;
    mxp1 = (m+2) / bpp;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    }
}  /* end of function load_int_el_split */


/*********************************************************************
 *  FUNCTION:  load_E_el_split -- loads the E elements of a collmat  *
 *  M = collmat                                                      *
 *  f = 0, use xl_yg numbering                                       *
 *  f = 1, use yl_xg numbering                                       *
 *  n, m = indices for which E element                               *
 *********************************************************************/
void load_E_el_split (M, fp, f, n, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, m, nproc, bpp, *ctr, size;
{
int i, j, k, nx, nxp1, mx, nslaves;

if (f == 0)  /* xl_yg */
    {
    nx   = (n+1) / bpp;
    nxp1 = (n+2) / bpp;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
	fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
        }
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nxp1+(++(ctr[nxp1]))]);
    }
else  /* f == 1, yl_xg */
    {
    mx = (m+1) / bpp;
    nslaves = nproc - 1;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    }
}  /* end of function load_E_el_split */


/**********************************************************************
 *  FUNCTION:  load_NW_el_split -- loads the NW element of a collmat  *
 *  M = collmat                                                       *
 *  f = 0, use xl_yg numbering                                        *
 *  f = 1, use yl_xg numbering                                        *
 *  n = index for which E element                                     *
 **********************************************************************/
void load_NW_el_split (M, fp, f, n, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, nproc, bpp, *ctr, size;
{
int i, j, k, nx, nslaves;

if (f == 0)  /* xl_yg */
    {
    nx = (n+1) / bpp;
    nslaves = nproc - 1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    
    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
else  /* f == 1, yl_xg */
    {
    /*  W equations  */
    for (i = 0; i < 2; i++)
        {
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[++(ctr[0])]);
        }

    /*  SE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  NE equation  */
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
}  /* end of function load_NW_el_split */


/*********************************************************************
 *  FUNCTION:  load_N_el_split -- loads the N elements of a collmat  *
 *  M = collmat                                                      *
 *  f = 0, use xl_yg numbering                                       *
 *  f = 1, use yl_xg numbering                                       *
 *  n, m = indices for which N element                               *
 *********************************************************************/
void load_N_el_split (M, fp, f, n, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, m, nproc, bpp, *ctr, size;
{
int i, j, nx, mx, mxp1, nslaves;

if (f == 0)  /* xl_yg */
    {
    nx = (n+1) / bpp;
    nslaves = nproc - 1;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    
    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
	fscanf(fp, "%lf", &M[size*nx+(++(ctr[nx]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
else  /* f == 1, yl_xg */
    {
    mx   = (m+1) / bpp;
    mxp1 = (m+2) / bpp;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);

    /*  SE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);

    /*  NE equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*mx+(++(ctr[mx]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*mxp1+(++(ctr[mxp1]))]);
    }
}  /*  end of function load_N_el_split  */


/**********************************************************************
 *  FUNCTION:  load_NE_el_split -- loads the NE element of a collmat  *
 *  M = collmat                                                       *
 *  f = 0, use xl_yg numbering                                        *
 *  f = 1, use yl_xg numbering                                        *
 *  n, m = indices for which element                                  *
 **********************************************************************/
void load_NE_el_split (M, fp, f, n, m, nproc, bpp, ctr, size)
double *M;
FILE *fp;
int f, n, m, nproc, bpp, *ctr, size;
{
int i, j, k, ix, nslaves;

nslaves = nproc - 1;

if (f == 0)  /* xl_yg */
    {
    ix = (n+1) / bpp;

    /*  SW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NW equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);

    /*  SE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
	fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
        }
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NE equation  */
    for (i = 0; i <= 1; i++)
        {
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
        }
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    for (i = 0; i <= 1; i++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    fscanf(fp, "%lf", &M[++(ctr[0])]);
    }
else  /* f == 1, yl_xg */
    {
    ix = (m+1) / bpp;

    /*  SW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  NW equation  */
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*ix+(++(ctr[ix]))]);
    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
    for (j = 0; j < 2; j++)
        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);

    /*  E equations  */
    for (i = 0; i < 2; i++)
        {
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	for (j = 0; j < 2; j++)
	    for (k = 0; k < 2; k++)
	        fscanf(fp, "%lf", &M[size*nslaves+(++(ctr[nslaves]))]);
	fscanf(fp, "%lf", &M[++(ctr[0])]);
	for (j = 0; j < 2; j++)
	    fscanf(fp, "%lf", &M[++(ctr[0])]);
        }
    }
}  /* end of function load_NE_el_split */


/*************************************************************************************
 *  FUNCTION:  file2split_collmat --  read entries for a collmat from a file, place  *
 *         them in vectors where each vector will be sent to a particular processor  *
 *************************************************************************************/
void file2split_collmat (V, filename, nproc, bpp, xm, ym, num, loc, glob, start, size)
double *V;
char filename[80];
int nproc, bpp, xm, ym, num, loc, glob, size;
{
int i, j, m, n, locm1, locm2, globm2, globm3, *ctr;
FILE *fp;

locm2 = loc-2;
locm1 = loc-1;
globm3 = glob-3;
globm2 = glob-2;

ctr = ivector(nproc);
for (i = 0; i < nproc; i++)
    ctr[i] = start - 1;

fp = fopen(filename, "r");

if (xm < ym)
    {
    if (num == 0)  /* xl_yg */
        {	    
	load_SW_el_split (V, fp, num, ctr, size);
	for (m = 1; m <= locm2; m++)
	    load_S_el_split (V, fp, num, m, nproc, bpp, ctr, size);
	load_SE_el_split (V, fp, num, locm1, nproc, bpp, ctr, size);
	for (n = 0; n <= globm3; n++)
	    {
	    load_W_el_split (V, fp, num, n, nproc, bpp, ctr, size);
	    for (m = 1; m <= locm2; m++)
	        load_int_el_split (V, fp, num, n, m, nproc, bpp, ctr, size);
	    load_E_el_split (V, fp, num, n, locm1, nproc, bpp, ctr, size);
	    }
	load_NW_el_split (V, fp, num, globm2, nproc, bpp, ctr, size);
	for (m = 1; m <= locm2; m++)
	    load_N_el_split (V, fp, num, globm2, m, nproc, bpp, ctr, size);
	load_NE_el_split (V, fp, num, globm2, locm1, nproc, bpp, ctr, size);	
        }
    else /* num == 1, yl_xg */
        {
	load_SW_el_split (V, fp, num, ctr, size);
	for (m = 0; m <= globm3; m++)
	    load_S_el_split (V, fp, num, m, nproc, bpp, ctr, size);
	load_SE_el_split (V, fp, num, globm2, nproc, bpp, ctr, size);
	for (n = 1; n <= locm2; n++)
	    {
	    load_W_el_split (V, fp, num, n, nproc, bpp, ctr, size);
	    for (m = 0; m <= globm3; m++)
	        load_int_el_split (V, fp, num, n, m, nproc, bpp, ctr, size);
	    load_E_el_split (V, fp, num, n, globm2, nproc, bpp, ctr, size);
	    }
	load_NW_el_split (V, fp, num, locm1, nproc, bpp, ctr, size);
	for (m = 0; m <= globm3; m++)
	    load_N_el_split (V, fp, num, locm1, m, nproc, bpp, ctr, size);
	load_NE_el_split (V, fp, num, locm1, globm2, nproc, bpp, ctr, size);
        }
    }
else  /* xm >= ym */
    {
    if (num == 0)  /* xl_yg */
        {           
	load_SW_el_split (V, fp, num, ctr, size);
	for (n = 0; n <= globm3; n++)
	    load_W_el_split (V, fp, num, n, nproc, bpp, ctr, size);
	load_NW_el_split (V, fp, num, globm2, nproc, bpp, ctr, size);
	for (m = 1; m <= locm2; m++)
	    {
	    load_S_el_split (V, fp, num, m, nproc, bpp, ctr, size);
	    for (n = 0; n <= globm3; n++)
	        load_int_el_split (V, fp, num, n, m, nproc, bpp, ctr, size);
	    load_N_el_split (V, fp, num, globm2, m, nproc, bpp, ctr, size);
	    }
	load_SE_el_split (V, fp, num, locm1, nproc, bpp, ctr, size); 
	for (n = 0; n <= globm3; n++)
	    load_E_el_split (V, fp, num, n, locm1, nproc, bpp, ctr, size);
	load_NE_el_split (V, fp, num, globm2, locm1, nproc, bpp, ctr, size);	
        }
    else  /* num == 1, yl_xg */
        {
	load_SW_el_split (V, fp, num, ctr, size);
	for (n = 1; n <= locm2; n++)
	    load_W_el_split (V, fp, num, n, nproc, bpp, ctr, size);
	load_NW_el_split (V, fp, num, locm1, nproc, bpp, ctr, size);
	for (m = 0; m <= globm3; m++)
	    {
	    load_S_el_split (V, fp, num, m, nproc, bpp, ctr, size);
	    for (n = 1; n <= locm2; n++)
	        load_int_el_split (V, fp, num, n, m, nproc, bpp, ctr, size);
	    load_N_el_split (V, fp, num, locm1, m, nproc, bpp, ctr, size);
	    }
	load_SE_el_split (V, fp, num, globm2, nproc, bpp, ctr, size);
	for (n = 1; n <= locm2; n++)
	    load_E_el_split (V, fp, num, n, globm2, nproc, bpp, ctr, size);
	load_NE_el_split (V, fp, num, locm1, globm2, nproc, bpp, ctr, size);
        }
    }

fclose(fp);

free(ctr);
}  /* end of function file2split_collmat */


/*********************************************************************
 *  FUNCTION:  unload_SW_M_n0 -- from the Master vector, unload the  *
 *             SW element using xl_yg numbering                      *
 *********************************************************************/
void unload_SW_M_n0 (M, V, ctr)
struct first_proc_collmat *M;
double *V;
int *ctr;
{
int i;

/*  SW equation  */
M->AF[0].a[0].el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->BF[0].a[0].el[0].el[0][i] = V[++(*ctr)];
M->AF[0].a[0].el[0][1] = V[++(*ctr)];
M->AF[0].b[0].el[0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].a[0].el[1].el[0][i] = V[++(*ctr)];
    M->BF[0].b[0]      .el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
M->CF[0].a[0].el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[0].a[0].A.el[0][i] = V[++(*ctr)];
M->CF[0].a[0].el[0][1] = V[++(*ctr)];
M->CF[0].b[0].el[0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].a[0].B    .el[0][i] = V[++(*ctr)];
    M->A[0].b[0].el[0].el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
M->AF[0].a[0].el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->BF[0].a[0].el[0].el[1][i] = V[++(*ctr)];
M->AF[0].a[0].el[1][1] = V[++(*ctr)];
M->AF[0].b[0].el[1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].a[0].el[1].el[1][i] = V[++(*ctr)];
    M->BF[0].b[0]      .el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
M->CF[0].a[0].el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[0].a[0].C.el[0][i] = V[++(*ctr)];
M->CF[0].a[0].el[1][1] = V[++(*ctr)];
M->CF[0].b[0].el[1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].a[0].D    .el[0][i] = V[++(*ctr)];
    M->A[0].b[0].el[1].el[0][i] = V[++(*ctr)];
    }
}  /* end of function unload_SW_M_n0 */


/******************************************************************
 *  FUNCTION:  unload_W_M_n0 -- from the Master vector, unload a  *
 *             W element using xl_yg numbering                    *
 ******************************************************************/
void unload_W_M_n0 (M, V, ctr, n, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, n, bppm2;
{
int i, np1 = n + 1;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    M->A[n].a[0].A.el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->B[n].a[0].el[0].el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[0].B    .el[1][i] = V[++(*ctr)];
    M->A[n].b[0].el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].a[0].el[1].el[0][i] = V[++(*ctr)];
    M->B[n].b[0]      .el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
for (i = 0; i <= 1; i++)
    M->C[n].a[0].el[0].el[0][i] = V[++(*ctr)];
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[0].A.el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->C[n].a[0].el[1].el[0][i] = V[++(*ctr)];
    M->C[n].b[0]      .el[0][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[0].B    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[0].el[0].el[0][i] = V[++(*ctr)];
        }

/*  SE equation  */
for (i = 0; i <= 1; i++)
    M->A[n].a[0].C.el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->B[n].a[0].el[0].el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[0].D    .el[1][i] = V[++(*ctr)];
    M->A[n].b[0].el[1].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].a[0].el[1].el[1][i] = V[++(*ctr)];
    M->B[n].b[0]      .el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
for (i = 0; i <= 1; i++)
    M->C[n].a[0].el[0].el[1][i] = V[++(*ctr)];
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[0].C.el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->C[n].a[0].el[1].el[1][i] = V[++(*ctr)];
    M->C[n].b[0]      .el[1][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[0].D    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[0].el[1].el[0][i] = V[++(*ctr)];
        }
}  /* end of function unload_W_M_n0 */


/*********************************************************************
 *  FUNCTION:  unload_NW_M_n0 -- from the Master vector, unload the  *
 *             NW element using xl_yg numbering                      *
 *********************************************************************/
void unload_NW_M_n0 (M, V, ctr)
struct first_proc_collmat *M;
double *V;
int *ctr;
{
/*  NW equation  */
M->AL[0].a[0].el[0][0] = V[++(*ctr)];
M->AL[0].a[0].el[0][1] = V[++(*ctr)];
M->AL[0].b[0].el[0] = V[++(*ctr)];

/*  NE equation  */
M->AL[0].a[0].el[1][0] = V[++(*ctr)];
M->AL[0].a[0].el[1][1] = V[++(*ctr)];
M->AL[0].b[0].el[1] = V[++(*ctr)];
}  /* end of function unload_NW_M_n0 */


/******************************************************************
 *  FUNCTION:  unload_S_M_n0 -- from the Master vector, unload a  *
 *             S element using xl_yg numbering                    *
 ******************************************************************/
void unload_S_M_n0 (M, V, ctr, m)
struct first_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, mm1;

mm1 = m - 1;
  
/*  SW equation  */
M->AF[0].c[mm1].el[0] = V[++(*ctr)];
M->AF[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->BF[0].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
M->AF[0].a[m].el[0][1] = V[++(*ctr)];
M->AF[0].b[m].el[0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].a[m].el[1].el[0][i] = V[++(*ctr)];
    M->BF[0].b[m]      .el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
M->CF[0].c[mm1].el[0] = V[++(*ctr)];
M->CF[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].c[mm1].el[0].el[0][i] = V[++(*ctr)];
    M->A[0].a[m].A.el[0][i] = V[++(*ctr)];
    }
M->CF[0].a[m].el[0][1] = V[++(*ctr)];
M->CF[0].b[m].el[0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].a[m].B    .el[0][i] = V[++(*ctr)];
    M->A[0].b[m].el[0].el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
M->AF[0].c[mm1].el[1] = V[++(*ctr)];
M->AF[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->BF[0].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
M->AF[0].a[m].el[1][1] = V[++(*ctr)];
M->AF[0].b[m].el[1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].a[m].el[1].el[1][i] = V[++(*ctr)];
    M->BF[0].b[m]      .el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
M->CF[0].c[mm1].el[1] = V[++(*ctr)];
M->CF[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].c[mm1].el[1].el[0][i] = V[++(*ctr)];
    M->A[0].a[m].C.el[0][i] = V[++(*ctr)];
    }
M->CF[0].a[m].el[1][1] = V[++(*ctr)];
M->CF[0].b[m].el[1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].a[m].D    .el[0][i] = V[++(*ctr)];
    M->A[0].b[m].el[1].el[0][i] = V[++(*ctr)];
    }
}  /* end of function unload_S_M_n0 */


/*********************************************************************
 *  FUNCTION:  unload_int_M_n0 -- from the Master vector, unload an  *
 *             interior element using xl_yg numbering                *
 *********************************************************************/
void unload_int_M_n0 (M, V, ctr, m, n, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, m, n, bppm2;
{  
int i, mm1, np1;

mm1 = m - 1;
np1 = n + 1;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
    M->A[n].b[m].el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];
    M->B[n].b[m]      .el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
for (i = 0; i <= 1; i++)
    {
    M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
    M->C[n].b[m]      .el[0][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[m].el[0].el[0][i] = V[++(*ctr)];
        }

/*  SE equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
    M->A[n].b[m].el[1].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];
    M->B[n].b[m]      .el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
for (i = 0; i <= 1; i++)
    {
    M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
    M->C[n].b[m]      .el[1][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[m].D    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[m].el[1].el[0][i] = V[++(*ctr)];
        }		
}  /* end of function unload_int_M_n0 */


/******************************************************************
 *  FUNCTION:  unload_N_M_n0 -- from the Master vector, unload a  *
 *             N element using xl_yg numbering                    *
 ******************************************************************/
void unload_N_M_n0 (M, V, ctr, m)
struct first_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, mm1;

mm1 = m - 1;
  
/*  NW equation  */
M->AL[0].c[mm1].el[0] = V[++(*ctr)];
M->AL[0].a[m].el[0][0] = V[++(*ctr)];
M->AL[0].a[m].el[0][1] = V[++(*ctr)];
M->AL[0].b[m].el[0] = V[++(*ctr)];

/*  NE equation  */
M->AL[0].c[mm1].el[1] = V[++(*ctr)];
M->AL[0].a[m].el[1][0] = V[++(*ctr)];
M->AL[0].a[m].el[1][1] = V[++(*ctr)];
M->AL[0].b[m].el[1] = V[++(*ctr)];
}  /* end of function unload_N_M_n0 */


/*******************************************************************
 *  FUNCTION:  unload_SE_M_n0 -- from the Master vector, unload a  *
 *             SE element using xl_yg numbering                    *
 *******************************************************************/
void unload_SE_M_n0 (M, V, ctr, m)
struct first_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, mm1;

mm1 = m - 1;

/*  SW equation  */
M->AF[0].c[mm1].el[0] = V[++(*ctr)];
M->AF[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->BF[0].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    } 
M->AF[0].a[m].el[0][1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->BF[0].a[m].el[1].el[0][i] = V[++(*ctr)];

/*  NW equation  */
M->CF[0].c[mm1].el[0] = V[++(*ctr)];
M->CF[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].c[mm1].el[0].el[0][i] = V[++(*ctr)];
    M->A[0].a[m].A.el[0][i] = V[++(*ctr)];
    } 
M->CF[0].a[m].el[0][1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[0].a[m].B    .el[0][i] = V[++(*ctr)];

/*  SE equation  */
M->AF[0].c[mm1].el[1] = V[++(*ctr)];
M->AF[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->BF[0].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->BF[0].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
M->AF[0].a[m].el[1][1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->BF[0].a[m].el[1].el[1][i] = V[++(*ctr)];

/*  NE equation  */
M->CF[0].c[mm1].el[1] = V[++(*ctr)];
M->CF[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[0].c[mm1].el[1].el[0][i] = V[++(*ctr)];
    M->A[0].a[m].C.el[0][i] = V[++(*ctr)];
    }
M->CF[0].a[m].el[1][1] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[0].a[m].D.el[0][i] = V[++(*ctr)];
}  /* end of function unload_SE_M_n0 */


/*******************************************************************
 *  FUNCTION:  unload_E_M_n0 -- from the Master vector, unload an  *
 *             E element using xl_yg numbering                     *
 *******************************************************************/
void unload_E_M_n0 (M, V, ctr, m, n, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, m, n, bppm2;
{
int i, mm1, np1;

mm1 = m - 1;
np1 = n + 1;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];

/*  NW equation  */
for (i = 0; i <= 1; i++)
    {
    M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];

/*  SE equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];

/*  NE equation  */
for (i = 0; i <= 1; i++)
    {
    M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
if (n < bppm2)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[m].D    .el[0][i] = V[++(*ctr)];
}  /* end of function unload_E_M_n0 */


/*******************************************************************
 *  FUNCTION:  unload_NE_M_n0 -- from the Master vector, unload a  *
 *             NE element using xl_yg numbering                    *
 *******************************************************************/
void unload_NE_M_n0 (M, V, ctr, m)
struct first_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, mm1;

mm1 = m - 1;

/*  NW equation  */
M->AL[0].c[mm1].el[0] = V[++(*ctr)];
M->AL[0].a[m]  .el[0][0] = V[++(*ctr)];
M->AL[0].a[m].el[0][1] = V[++(*ctr)];

/*  NE equation  */
M->AL[0].c[mm1].el[1] = V[++(*ctr)];
M->AL[0].a[m]  .el[1][0] = V[++(*ctr)];
M->AL[0].a[m].el[1][1] = V[++(*ctr)];
}  /* end of function unload_NE_M_n0 */


/*********************************************************************
 *  FUNCTION:  unload_SW_M_n1 -- from the Master vector, unload the  *
 *             SW element using yl_xg numbering                      *
 *********************************************************************/
void unload_SW_M_n1 (M, V, ctr)
struct first_proc_collmat *M;
double *V;
int *ctr;
{
int i, j, k;

/*  W equations  */
for (i = 0; i < 2; i++)
    {
    for (j = 0; j < 2; j++)
        M->AF[0].a[0].el[i][j] = V[++(*ctr)];
    M->AF[0].b[0].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->BF[0].a[0].el[j].el[i][k] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->BF[0].b[0].el[i][j] = V[++(*ctr)];
    }

/*  SE equation  */
for (j = 0; j < 2; j++)
    M->CF[0].a[0].el[0][j] = V[++(*ctr)];
M->CF[0].b[0].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[0].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[0].B.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].b[0].el[0].el[0][j] = V[++(*ctr)];

/*  NE equation  */
for (j = 0; j < 2; j++)
    M->CF[0].a[0].el[1][j] = V[++(*ctr)];
M->CF[0].b[0].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[0].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[0].D.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].b[0].el[1].el[0][j] = V[++(*ctr)];
}  /* end of function unload_SW_M_n1 */


/******************************************************************
 *  FUNCTION:  unload_W_M_n1 -- from the Master vector, unload a  *
 *             W element using yl_xg numbering                    *
 ******************************************************************/
void unload_W_M_n1 (M, V, ctr, n)
struct first_proc_collmat *M;
double *V;
int *ctr, n;
{
int i, j, k, nm1;

nm1 = n - 1;

/*  W equations  */
for (i = 0; i < 2; i++)
    {
    M->AF[0].c[nm1].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->AF[0].a[n].el[i][j] = V[++(*ctr)];
    M->AF[0].b[n].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->BF[0].c[nm1].el[i][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->BF[0].a[n].el[j].el[i][k] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->BF[0].b[n].el[i][j] = V[++(*ctr)];
    }

/*  SE equation  */
M->CF[0].c[nm1].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->CF[0].a[n].el[0][j] = V[++(*ctr)];
M->CF[0].b[n].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].c[nm1].el[0].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].B.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].b[n].el[0].el[0][j] = V[++(*ctr)];

/*  NE equation  */
M->CF[0].c[nm1].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->CF[0].a[n].el[1][j] = V[++(*ctr)];
M->CF[0].b[n].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].c[nm1].el[1].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].D.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].b[n].el[1].el[0][j] = V[++(*ctr)];
}  /* end of function unload_W_M_n1 */


/*******************************************************************
 *  FUNCTION:  unload_NW_M_n1 -- from the Master vector, unload a  *
 *             NW element using yl_xg numbering                    *
 *******************************************************************/
void unload_NW_M_n1 (M, V, ctr, n)
struct first_proc_collmat *M;
double *V;
int *ctr, n;
{
int i, j, k, nm1;

nm1 = n - 1;

/*  W equations  */
for (i = 0; i < 2; i++)
    {
    M->AF[0].c[nm1].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->AF[0].a[n].el[i][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->BF[0].c[nm1].el[i][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->BF[0].a[n].el[j].el[i][k] = V[++(*ctr)];
    }

/*  SE equation  */
M->CF[0].c[nm1].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->CF[0].a[n].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].c[nm1].el[0].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].B.el[0][j] = V[++(*ctr)];

/*  NE equation  */
M->CF[0].c[nm1].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->CF[0].a[n].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].c[nm1].el[1].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[0].a[n].D.el[0][j] = V[++(*ctr)];
}  /* end of function unload_NW_M_n1 */


/*******************************************************************
 *  FUNCTION:  unload_S_M_n1 -- from the Master vector, unload a   *
 *             S element using yl_xg numbering                     *
 *******************************************************************/
void unload_S_M_n1 (M, V, ctr, m, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, m, bppm2;
{
int i, j, mp1;

mp1 = m + 1;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].a[0].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[0].B.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[0].el[0].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[0].el[i].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].b[0].el[0][j] = V[++(*ctr)];

/*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].a[0].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[0].D.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[0].el[1].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[0].el[i].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].b[0].el[1][j] = V[++(*ctr)];

/*  SE equation  */
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[0].el[i].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->C[m].b[0].el[0][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].A.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].B.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[0].el[0].el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[0].el[i].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->C[m].b[0].el[1][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].C.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].D.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[0].el[1].el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_S_M_n1 */


/*********************************************************************
 *  FUNCTION:  unload_int_M_n1 -- from the Master vector, unload an  *
 *             interior element using yl_xg numbering                *
 *********************************************************************/
unload_int_M_n1 (M, V, ctr, n, m, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, n, m, bppm2;
{
int i, j, mp1, nm1;

mp1 = m + 1;
nm1 = n - 1;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[n].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].b[n].el[0][j] = V[++(*ctr)];

/*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[n].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].b[n].el[1][j] = V[++(*ctr)];

/*  SE equation  */
for (j = 0; j < 2; j++)
    M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->C[m].b[n].el[0][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[n].el[0].el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
for (j = 0; j < 2; j++)
    M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->C[m].b[n].el[1][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[n].el[1].el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_int_M_n1 */


/*******************************************************************
 *  FUNCTION:  unload_N_M_n1 -- from the Master vector, unload an  *
 *             N element using yl_xg numbering                     *
 *******************************************************************/
unload_N_M_n1 (M, V, ctr, n, m, bppm2)
struct first_proc_collmat *M;
double *V;
int *ctr, n, m, bppm2;
{
int i, j, mp1, nm1;

mp1 = m + 1;
nm1 = n - 1;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];

/*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];

/*  SE equation  */
for (j = 0; j < 2; j++)
    M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
for (j = 0; j < 2; j++)
    M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
        M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
if (m < bppm2)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_N_M_n1 */


/********************************************************************
 *  FUNCTION:  unload_SE_M_n1 -- from the Master vector, unload an  *
 *             SE element using yl_xg numbering                     *
 ********************************************************************/
unload_SE_M_n1 (M, V, ctr)
struct first_proc_collmat *M;
double *V;
int *ctr;
{
int i, j;

for (i = 0; i < 2; i++)
    {
    for (j = 0; j < 2; j++)
        M->AL[0].a[0].el[i][j] = V[++(*ctr)];
    M->AL[0].b[0].el[i] = V[++(*ctr)];
    }
}  /* end of function unload_SE_M_n1 */


/*******************************************************************
 *  FUNCTION:  unload_E_M_n1 -- from the Master vector, unload an  *
 *             E element using yl_xg numbering                     *
 *******************************************************************/
unload_E_M_n1 (M, V, ctr, n)
struct first_proc_collmat *M;
double *V;
int *ctr, n;
{
int i, j, nm1;

nm1 = n - 1;

for (i = 0; i < 2; i++)
    {
    M->AL[0].c[nm1].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->AL[0].a[n].el[i][j] = V[++(*ctr)];
    M->AL[0].b[n].el[i] = V[++(*ctr)];
    }
}  /* end of function unload_E_M_n1 */


/********************************************************************
 *  FUNCTION:  unload_NE_M_n1 -- from the Master vector, unload an  *
 *             NE element using yl_xg numbering                     *
 ********************************************************************/
unload_NE_M_n1 (M, V, ctr, n)
struct first_proc_collmat *M;
double *V;
int *ctr, n;
{
int i, j, nm1;

nm1 = n - 1;

for (i = 0; i < 2; i++)
    {
    M->AL[0].c[nm1].el[i] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->AL[0].a[n].el[i][j] = V[++(*ctr)];
    }
}  /* end of function unload_NE_M_n1 */


/*******************************
 *  FUNCTION:  load_collmat_M  *
 *******************************/
void load_collmat_M (M, V, bpp, xm, ym, num, loc, start)
struct first_proc_collmat *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int ctr, n, m, np1, mm1, locm2, locm1, i, bppm2;

locm2  = loc  - 2;
locm1  = loc  - 1;
bppm2  = bpp  - 2;
ctr = start - 1;

if (xm < ym)
    {
    if (num == 0)  /* xl_yg */
        {	    
	unload_SW_M_n0 (M, V, &ctr);
	for (m = 1; m <= locm2; m++)
	    unload_S_M_n0 (M, V, &ctr, m);
	unload_SE_M_n0 (M, V, &ctr, locm1);
	for (n = 0; n <= bppm2; n++)
	    {
	    unload_W_M_n0 (M, V, &ctr, n, bppm2);
	    for (m = 1; m <= locm2; m++)
	        unload_int_M_n0 (M, V, &ctr, m, n, bppm2);
	    unload_E_M_n0 (M, V, &ctr, locm1, n, bppm2);
	    }
	unload_NW_M_n0 (M, V, &ctr);
	for (m = 1; m <= locm2; m++)
	    unload_N_M_n0 (M, V, &ctr, m);
	unload_NE_M_n0 (M, V, &ctr, locm1);
        }
    else /* num == 1, yl_xg */
        {
	unload_SW_M_n1 (M, V, &ctr);
	for (m = 0; m <= bppm2; m++)
	    unload_S_M_n1 (M, V, &ctr, m, bppm2);
	unload_SE_M_n1 (M, V, &ctr);
	for (n = 1; n <= locm2; n++)
	    {
	    unload_W_M_n1 (M, V, &ctr, n);
	    for (m = 0; m <= bppm2; m++)
	        unload_int_M_n1 (M, V, &ctr, n, m, bppm2);
	    unload_E_M_n1 (M, V, &ctr, n);
	    }
	unload_NW_M_n1 (M, V, &ctr, locm1);
	for (m = 0; m <= bppm2; m++)
	    unload_N_M_n1 (M, V, &ctr, locm1, m, bppm2);
	unload_NE_M_n1 (M, V, &ctr, locm1);
        }
    }
else  /* xm >= ym */
    {
    if (num == 0)  /* xl_yg */
        {    
	unload_SW_M_n0 (M, V, &ctr);
	for (n = 0; n <= bppm2; n++)
	    unload_W_M_n0 (M, V, &ctr, n, bppm2);
	unload_NW_M_n0 (M, V, &ctr);
	for (m = 1; m <= locm2; m++)
	    {
	    unload_S_M_n0 (M, V, &ctr, m);
	    for (n = 0; n <= bppm2; n++)
		unload_int_M_n0 (M, V, &ctr, m, n, bppm2);
	    unload_N_M_n0 (M, V, &ctr, m);
	    }
	m = locm1;
	unload_SE_M_n0 (M, V, &ctr, m);
	for (n = 0; n <= bppm2; n++)
	    unload_E_M_n0 (M, V, &ctr, m, n, bppm2);
	unload_NE_M_n0 (M, V, &ctr, m);
        }
    else  /* num == 1, yl_xg */
        {
	unload_SW_M_n1 (M, V, &ctr);
	for (n = 1; n <= locm2; n++)
	    unload_W_M_n1 (M, V, &ctr, n);
	unload_NW_M_n1 (M, V, &ctr, locm1);
	for (m = 0; m <= bppm2; m++)
	    {
	    unload_S_M_n1 (M, V, &ctr, m, bppm2);
	    for (n = 1; n <= locm2; n++)
	        unload_int_M_n1 (M, V, &ctr, n, m, bppm2);
	    unload_N_M_n1 (M, V, &ctr, locm1, m, bppm2);
	    }
	unload_SE_M_n1 (M, V, &ctr);
	for (n = 1; n <= locm2; n++)
	    unload_E_M_n1 (M, V, &ctr, n);
	unload_NE_M_n1 (M, V, &ctr, locm1);
        }
    }

}  /* end of function load_collmat_M */


/****************************************************************
 *  FUNCTION:  unload_W_S_n0 -- from a slave vector, unload an  *
 *             W element using xl_yg numbering                  *
 ****************************************************************/
unload_W_S_n0 (M, V, ctr, n, bpp)
struct proc_collmat *M;
double *V;
int *ctr, n, bpp;
{
int i, bppm1, np1;

bppm1 = bpp - 1;
np1 = n + 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        M->A[n].a[0].A.el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[0].el[0].el[0][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[0].B    .el[1][i] = V[++(*ctr)];
	M->A[n].b[0].el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[0].el[1].el[0][i] = V[++(*ctr)];
	M->B[n].b[0]      .el[0][i] = V[++(*ctr)];
        }
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[0].el[0].el[0][i] = V[++(*ctr)];
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[0].A.el[0][i] = V[++(*ctr)];
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[0].el[1].el[0][i] = V[++(*ctr)];
	M->C[n].b[0]      .el[0][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[0].B    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[0].el[0].el[0][i] = V[++(*ctr)];
        }

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        M->A[n].a[0].C.el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[0].el[0].el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[0].D    .el[1][i] = V[++(*ctr)];
	M->A[n].b[0].el[1].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[0].el[1].el[1][i] = V[++(*ctr)];
	M->B[n].b[0]      .el[1][i] = V[++(*ctr)];
        }
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[0].el[0].el[1][i] = V[++(*ctr)];
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[0].C.el[0][i] = V[++(*ctr)];
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[0].el[1].el[1][i] = V[++(*ctr)];
	M->C[n].b[0]      .el[1][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[0].D    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[0].el[1].el[0][i] = V[++(*ctr)];
        }
}  /* end of function unload_W_S_n0 */  


/******************************************************************
 *  FUNCTION:  unload_int_S_n0 -- from a slave vector, unload an  *
 *             interior element using xl_yg numbering             *
 ******************************************************************/
unload_int_S_n0 (M, V, ctr, m, n, bpp)
struct proc_collmat *M;
double *V;
int *ctr, m, n, bpp;
{
int i, np1, mm1, bppm1;

np1 = n+1;
mm1 = m-1;
bppm1 = bpp - 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
	M->A[n].b[m].el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];
	M->B[n].b[m]      .el[0][i] = V[++(*ctr)];
        }
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
        }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
	M->C[n].b[m]      .el[0][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[m].el[0].el[0][i] = V[++(*ctr)];
        }

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
	M->A[n].b[m].el[1].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];
	M->B[n].b[m]      .el[1][i] = V[++(*ctr)];
        }
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
        }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
	M->C[n].b[m]      .el[1][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].a[m].D    .el[0][i] = V[++(*ctr)];
	M->A[np1].b[m].el[1].el[0][i] = V[++(*ctr)];
        }
}  /* end of function unload_int_S_n0 */


/****************************************************************
 *  FUNCTION:  unload_E_S_n0 -- from a slave vector, unload an  *
 *             E element using xl_yg numbering                  *
 ****************************************************************/
unload_E_S_n0 (M, V, ctr, m, n, bpp)
struct proc_collmat *M;
double *V;
int *ctr, m, n, bpp;
{
int i, mm1, np1, bppm1;

mm1 = m - 1;
np1 = n + 1;
bppm1 = bpp - 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
        }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        {
	M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
	M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
        }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
if (n < bppm1)
    for (i = 0; i <= 1; i++)
        M->A[np1].a[m].D  .el[0][i] = V[++(*ctr)];
}  /* end of function unload_E_S_n0 */


/****************************************************************
 *  FUNCTION:  unload_S_S_n1 -- from a slave vector, unload an  *
 *             S element using yl_xg numbering                  *
 ****************************************************************/
unload_S_S_n1 (M, V, ctr, m, bpp)
struct proc_collmat *M;
double *V;
int *ctr, m, bpp;
{
int i, j, mp1, bppm1;

mp1 = m + 1;
bppm1 = bpp - 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].a[0].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[0].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[0].el[0].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[0].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[0].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].a[0].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[0].D.el[1][j] = V[++(*ctr)];
    for (j =0; j < 2; j++)
        M->A[m].b[0].el[1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[0].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[0].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[0].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[0].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].A.el[0][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].B.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[0].el[0].el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
if (m >= 0)
    {
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[0].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].C.el[0][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].a[0].D.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[0].el[1].el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_S_S_n1 */


/******************************************************************
 *  FUNCTION:  unload_int_S_n1 -- from a slave vector, unload an  *
 *             interior element using yl_xg numbering             *
 ******************************************************************/
unload_int_S_n1 (M, V, ctr, n, m, bpp)
struct proc_collmat *M;
double *V;
int *ctr, n, m, bpp;
{
int i, j, mp1, nm1, bppm1;

mp1 = m + 1;
nm1 = n - 1;
bppm1 = bpp - 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[n].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[n].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[n].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[n].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[n].el[0][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[n].el[0].el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[n].el[1][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].b[n].el[1].el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_int_S_n1 */


/****************************************************************
 *  FUNCTION:  unload_N_S_n1 -- from a slave vector, unload an  *
 *             N element using yl_xg numbering                  *
 ****************************************************************/
unload_N_S_n1 (M, V, ctr, n, m, bpp)
struct proc_collmat *M;
double *V;
int *ctr, n, m, bpp;
{
int i, j, nm1, mp1, bppm1;

nm1 = n - 1;
mp1 = m + 1;
bppm1 = bpp - 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
       M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];
    }

/*  NE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    }
if (m < bppm1)
    {
    for (j = 0; j < 2; j++)
        M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
    }
}  /* end of function unload_N_S_n1 */


/*******************************
 *  FUNCTION:  load_collmat_S  *
 *******************************/
void load_collmat_S (M, V, bpp, xm, ym, num, loc, start)
struct proc_collmat *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int ctr, n, m, np1, mm1, locm2, locm1, i, bppm2;

locm2  = loc  - 2;
locm1  = loc  - 1;
bppm2  = bpp  - 2;
ctr = start - 1;

if (xm < ym)
    {
    if (num == 0)  /* xl_yg */
	for (n = -1; n < bpp; n++)
	    {
	    unload_W_S_n0 (M, V, &ctr, n, bpp);
	    for (m = 1; m <= locm2; m++)
	        unload_int_S_n0 (M, V, &ctr, m, n, bpp);
	    unload_E_S_n0 (M, V, &ctr, locm1, n, bpp);
	    }
    else /* num == 1, yl_xg */
        {
	for (m = -1; m < bpp; m++)
	    unload_S_S_n1 (M, V, &ctr, m, bpp);
	for (n = 1; n <= locm2; n++)
	    for (m = -1; m < bpp; m++)
	        unload_int_S_n1 (M, V, &ctr, n, m, bpp);
	for (m = -1; m < bpp; m++)
	    unload_N_S_n1 (M, V, &ctr, locm1, m, bpp);
        }
    }
else  /* xm >= ym */
    {
    if (num == 0)  /* xl_yg */
        { 
	for (n = -1; n < bpp; n++)
	    unload_W_S_n0 (M, V, &ctr, n, bpp);
	for (m = 1; m <= locm2; m++)
	    for (n = -1; n < bpp; n++)
		unload_int_S_n0 (M, V, &ctr, m, n, bpp);
	for (n = -1; n < bpp; n++)
	    unload_E_S_n0 (M, V, &ctr, locm1, n, bpp);
        }
    else  /* num == 1, yl_xg */
        {
	for (m = -1; m < bpp; m++)
	    {
	    unload_S_S_n1 (M, V, &ctr, m, bpp);
	    for (n = 1; n <= locm2; n++)
	        unload_int_S_n1 (M, V, &ctr, n, m, bpp);
	    unload_N_S_n1 (M, V, &ctr, locm1, m, bpp);
	    }
        }
    }
}  /* end of function load_collmat_S */


/****************************************************************
 *  FUNCTION:  unload_W_L_n0 -- from a lslave vector, unload a  *
 *             W element using xl_yg numbering                  *
 ****************************************************************/
unload_W_L_n0 (M, V, ctr, n)
struct last_proc_collmat *M;
double *V;
int *ctr, n;
{
int i, np1;

np1 = n + 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        M->A[n].a[0].A.el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[0].el[0].el[0][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[0].B    .el[1][i] = V[++(*ctr)];
	M->A[n].b[0].el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[0].el[1].el[0][i] = V[++(*ctr)];
	M->B[n].b[0]      .el[0][i] = V[++(*ctr)];
        }
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[0].el[0].el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[np1].a[0].A.el[0][i] = V[++(*ctr)];
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[0].el[1].el[0][i] = V[++(*ctr)];
	M->C[n].b[0]      .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].a[0].B    .el[0][i] = V[++(*ctr)];
    M->A[np1].b[0].el[0].el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        M->A[n].a[0].C.el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[0].el[0].el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[0].D    .el[1][i] = V[++(*ctr)];
	M->A[n].b[0].el[1].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[0].el[1].el[1][i] = V[++(*ctr)];
	M->B[n].b[0]      .el[1][i] = V[++(*ctr)];
        }
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[0].el[0].el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[np1].a[0].C.el[0][i] = V[++(*ctr)];
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[0].el[1].el[1][i] = V[++(*ctr)];
	M->C[n].b[0]      .el[1][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].a[0].D    .el[0][i] = V[++(*ctr)];
    M->A[np1].b[0].el[1].el[0][i] = V[++(*ctr)];
    }
}  /* end of function unload_W_L_n0 */     


/*****************************************************************
 *  FUNCTION:  unload_NW_L_n0 -- from a lslave vector, unload a  *
 *             NW element using xl_yg numbering                  *
 *****************************************************************/
unload_NW_L_n0 (M, V, ctr, n)
struct last_proc_collmat *M;
double *V;
int *ctr, n;
{
int i;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    M->A[n].a[0].A.el[1][i] = V[++(*ctr)];
M->BL[0].a[0].el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[0].B    .el[1][i] = V[++(*ctr)];
    M->A[n].b[0].el[0].el[1][i] = V[++(*ctr)];
    }
M->BL[0].a[0].el[0][1] = V[++(*ctr)];
M->BL[0].b[0].el[0] = V[++(*ctr)];

/*  NW equation  */
for (i = 0; i <= 1; i++)
    M->CL[0].a[0].el[0].el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->CL[0].a[0].el[1].el[0][i] = V[++(*ctr)];
    M->CL[0].b[0]      .el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
for (i = 0; i <= 1; i++)
    M->A[n].a[0].C.el[1][i] = V[++(*ctr)];
M->BL[0].a[0].el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[0].D    .el[1][i] = V[++(*ctr)];
    M->A[n].b[0].el[1].el[1][i] = V[++(*ctr)];
    }
M->BL[0].a[0].el[1][1] = V[++(*ctr)];
M->BL[0].b[0].el[1] = V[++(*ctr)];

/*  NE equation  */
for (i = 0; i <= 1; i++)
    M->CL[0].a[0].el[0].el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->CL[0].a[0].el[1].el[1][i] = V[++(*ctr)];
    M->CL[0].b[0]      .el[1][i] = V[++(*ctr)];
    }
}  /* end of function unload_NW_L_n0 */


/******************************************************************
 *  FUNCTION:  unload_int_L_n0 -- from a lslave vector, unload a  *
 *             interior element using xl_yg numbering             *
 ******************************************************************/
unload_int_L_n0 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, mm1, np1;

mm1 = m - 1;
np1 = n + 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
	M->A[n].b[m].el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];
	M->B[n].b[m]      .el[0][i] = V[++(*ctr)];
        }
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
    M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
    }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
	M->C[n].b[m]      .el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];
    M->A[np1].b[m].el[0].el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
	M->A[n].b[m].el[1].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];
	M->B[n].b[m]      .el[1][i] = V[++(*ctr)];
        }
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
    M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
    }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
	M->C[n].b[m]      .el[1][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].a[m].D    .el[0][i] = V[++(*ctr)];
    M->A[np1].b[m].el[1].el[0][i] = V[++(*ctr)];
    }
}  /* end of function unload_int_L_n0 */


/******************************************************************
 *  FUNCTION:  unload_N_L_n0 -- from a lslave vector, unload a  *
 *             N element using xl_yg numbering                 *
 **************************************************************/
unload_N_L_n0 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, mm1;

mm1 = m - 1;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
    }
M->BL[0].c[mm1].el[0] = V[++(*ctr)];
M->BL[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
    M->A[n].b[m].el[0].el[1][i] = V[++(*ctr)];
    }
M->BL[0].a[m].el[0][1] = V[++(*ctr)];
M->BL[0].b[m].el[0] = V[++(*ctr)];

/*  NW equation  */
for (i = 0; i <= 1; i++)
    {
    M->CL[0].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->CL[0].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->CL[0].a[m].el[1].el[0][i] = V[++(*ctr)];
    M->CL[0].b[m]      .el[0][i] = V[++(*ctr)];
    }

/*  SE equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
    }
M->BL[0].c[mm1].el[1] = V[++(*ctr)];
M->BL[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    {
    M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
    M->A[n].b[m].el[1].el[1][i] = V[++(*ctr)];
    }
M->BL[0].a[m].el[1][1] = V[++(*ctr)];
M->BL[0].b[m].el[1] = V[++(*ctr)];

/*  NE equation  */
for (i = 0; i <= 1; i++)
    {
    M->CL[0].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->CL[0].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    {
    M->CL[0].a[m].el[1].el[1][i] = V[++(*ctr)];
    M->CL[0].b[m]      .el[1][i] = V[++(*ctr)];
    }
}  /* end of function unload_N_L_n0 */


/****************************************************************
 *  FUNCTION:  unload_E_L_n0 -- from a lslave vector, unload a  *
 *             E element using xl_yg numbering                  *
 ****************************************************************/
unload_E_L_n0 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, mm1, np1;

mm1 = m - 1;
np1 = n + 1;

/*  SW equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[m].el[1].el[0][i] = V[++(*ctr)];
    }

/*  NW equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[0][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[0][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].c[mm1].el[0].el[0][i] = V[++(*ctr)];
    M->A[np1].a[m]  .A    .el[0][i] = V[++(*ctr)];
    }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[m].el[1].el[0][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[np1].a[m].B    .el[0][i] = V[++(*ctr)];

/*  SE equation  */
if (n >= 0)
    {
    for (i = 0; i <= 1; i++)
        {
	M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
	M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        {
	M->B[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->B[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
    for (i = 0; i <= 1; i++)
        M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
    for (i = 0; i <= 1; i++)
        M->B[n].a[m].el[1].el[1][i] = V[++(*ctr)];
    }

/*  NE equation  */
if (n >= 0)
    for (i = 0; i <= 1; i++)
        {
	M->C[n].c[mm1]      .el[1][i] = V[++(*ctr)];
	M->C[n].a[m]  .el[0].el[1][i] = V[++(*ctr)];
        }
for (i = 0; i <= 1; i++)
    {
    M->A[np1].c[mm1].el[1].el[0][i] = V[++(*ctr)];
    M->A[np1].a[m]  .C    .el[0][i] = V[++(*ctr)];
    }
if (n >= 0)
    for (i = 0; i <= 1; i++)
        M->C[n].a[m].el[1].el[1][i] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[np1].a[m].D    .el[0][i] = V[++(*ctr)];
      
}  /* end of function unload_E_L_n0 */


/*****************************************************************
 *  FUNCTION:  unload_NE_L_n0 -- from a lslave vector, unload a  *
 *             NE element using xl_yg numbering                  *
 *****************************************************************/
unload_NE_L_n0 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, mm1;

mm1 = m - 1;

/*  SW equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[0].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .A    .el[1][i] = V[++(*ctr)];
    }
M->BL[0].c[mm1].el[0] = V[++(*ctr)];
M->BL[0].a[m]  .el[0][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[n].a[m].B    .el[1][i] = V[++(*ctr)];
M->BL[0].a[m].el[0][1] = V[++(*ctr)];

/*  NW equation  */
for (i = 0; i <= 1; i++)
    {
    M->CL[0].c[mm1]      .el[0][i] = V[++(*ctr)];
    M->CL[0].a[m]  .el[0].el[0][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    M->CL[0].a[m].el[1].el[0][i] = V[++(*ctr)];

/*  SE equation  */
for (i = 0; i <= 1; i++)
    {
    M->A[n].c[mm1].el[1].el[1][i] = V[++(*ctr)];
    M->A[n].a[m]  .C    .el[1][i] = V[++(*ctr)];
    }
M->BL[0].c[mm1].el[1] = V[++(*ctr)];
M->BL[0].a[m]  .el[1][0] = V[++(*ctr)];
for (i = 0; i <= 1; i++)
    M->A[n].a[m].D    .el[1][i] = V[++(*ctr)];
M->BL[0].a[m].el[1][1] = V[++(*ctr)];

/*  NE equation  */
for (i = 0; i <= 1; i++)
    {
    M->CL[0].c[mm1]      .el[1][i] = V[++(*ctr)];
    M->CL[0].a[m]  .el[0].el[1][i] = V[++(*ctr)];
    }
for (i = 0; i <= 1; i++)
    M->CL[0].a[m].el[1].el[1][i] = V[++(*ctr)];
}  /* end of function unload_NE_L_n0 */


/****************************************************************
 *  FUNCTION:  unload_S_L_n1 -- from a lslave vector, unload a  *
 *             S element using yl_xg numbering                  *
 ****************************************************************/
unload_S_L_n1 (M, V, ctr, m)
struct last_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, j, mp1;

mp1 = m + 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].a[0].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[0].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[0].el[0].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[0].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[0].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].a[0].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[0].D.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[0].el[1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[0].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[0].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[0].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[0].el[0][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)
    M->A[mp1].a[0].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[0].B.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].b[0].el[0].el[0][j] = V[++(*ctr)];

/*  NE equation  */
if (m >= 0)
    {
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[0].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[0].el[1][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)	
    M->A[mp1].a[0].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[0].D.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].b[0].el[1].el[0][j] = V[++(*ctr)];
}  /* end of function unload_S_L_n1 */


/******************************************************************
 *  FUNCTION:  unload_int_L_n1 -- from a lslave vector, unload a  *
 *             interior element using yl_xg numbering             *
 ******************************************************************/
unload_int_L_n1 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, j, nm1, mp1;

nm1 = n - 1;
mp1 = m + 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[n].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[n].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].b[n].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].b[n].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[n].el[0][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)
    M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].b[n].el[0].el[0][j] = V[++(*ctr)];

/*  NE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->C[m].b[n].el[1][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)
    M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].b[n].el[1].el[0][j] = V[++(*ctr)];
}  /* end of function unload_int_L_n1 */


/****************************************************************
 *  FUNCTION:  unload_N_L_n1 -- from a lslave vector, unload a  *
 *             N element using yl_xg numbering                  *
 ****************************************************************/
unload_N_L_n1 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, j, mp1, nm1;

mp1 = m + 1;
nm1 = n - 1;

/*  SW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    }

/*  NW equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->B[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->B[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    }

/*  SE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[0][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[0][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)
    M->A[mp1].c[nm1].el[0].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].A.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].B.el[0][j] = V[++(*ctr)];

/*  NE equation  */
if (m >= 0)
    {
    for (j = 0; j < 2; j++)
        M->C[m].c[nm1].el[1][j] = V[++(*ctr)];
    for (i = 0; i < 2; i++)
        for (j = 0; j < 2; j++)
	    M->C[m].a[n].el[i].el[1][j] = V[++(*ctr)];
    }
for (j = 0; j < 2; j++)
    M->A[mp1].c[nm1].el[1].el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].C.el[0][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[mp1].a[n].D.el[0][j] = V[++(*ctr)];
}  /* end of function unload_N_L_n1 */


/*****************************************************************
 *  FUNCTION:  unload_SE_L_n1 -- from a lslave vector, unload a  *
 *             SE element using yl_xg numbering                  *
 *****************************************************************/
unload_SE_L_n1 (M, V, ctr, m)
struct last_proc_collmat *M;
double *V;
int *ctr, m;
{
int i, j, k;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].a[0].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[0].B.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
     M->A[m].b[0].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[0].el[0][j] = V[++(*ctr)];
M->BL[0].b[0].el[0] = V[++(*ctr)];

/*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].a[0].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[0].D.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[0].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[0].el[1][j] = V[++(*ctr)];
M->BL[0].b[0].el[1] = V[++(*ctr)];

/*  E equations  */
for (i = 0; i < 2; i++)
    {
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->CL[0].a[0].el[j].el[i][k] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->CL[0].b[0].el[i][j] = V[++(*ctr)];
    }
}  /* end of function unload_SE_L_n1 */


/****************************************************************
 *  FUNCTION:  unload_E_L_n1 -- from a lslave vector, unload a  *
 *             E element using yl_xg numbering                  *
 ****************************************************************/
unload_E_L_n1 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, j, k, nm1;

nm1 = n - 1;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[n].el[0].el[1][j] = V[++(*ctr)];
M->BL[0].c[nm1].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[n].el[0][j] = V[++(*ctr)];
M->BL[0].b[n].el[0] = V[++(*ctr)];

          /*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].b[n].el[1].el[1][j] = V[++(*ctr)];
M->BL[0].c[nm1].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[n].el[1][j] = V[++(*ctr)];
M->BL[0].b[n].el[1] = V[++(*ctr)];

/*  E equations  */
for (i = 0; i < 2; i++)
    {
    for (j = 0; j < 2; j++)
        M->CL[0].c[nm1].el[i][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->CL[0].a[n].el[j].el[i][k] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        M->CL[0].b[n].el[i][j] = V[++(*ctr)];
    }
}  /* end of function unload_E_L_n1 */


/*****************************************************************
 *  FUNCTION:  unload_NE_L_n1 -- from a lslave vector, unload a  *
 *             NE element using yl_xg numbering                  *
 *****************************************************************/
unload_NE_L_n1 (M, V, ctr, n, m)
struct last_proc_collmat *M;
double *V;
int *ctr, n, m;
{
int i, j, k, nm1;

nm1 = n - 1;

/*  SW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[0].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].A.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].B.el[1][j] = V[++(*ctr)];
M->BL[0].c[nm1].el[0] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[n].el[0][j] = V[++(*ctr)];

/*  NW equation  */
for (j = 0; j < 2; j++)
    M->A[m].c[nm1].el[1].el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].C.el[1][j] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->A[m].a[n].D.el[1][j] = V[++(*ctr)];
M->BL[0].c[nm1].el[1] = V[++(*ctr)];
for (j = 0; j < 2; j++)
    M->BL[0].a[n].el[1][j] = V[++(*ctr)];

/*  E equations  */
for (i = 0; i < 2; i++)
    {
    for (j = 0; j < 2; j++)
        M->CL[0].c[nm1].el[i][j] = V[++(*ctr)];
    for (j = 0; j < 2; j++)
        for (k = 0; k < 2; k++)
	    M->CL[0].a[n].el[j].el[i][k] = V[++(*ctr)];
    }
}  /* end of function unload_NE_L_n1 */


/*******************************
 *  FUNCTION:  load_collmat_L  *
 *******************************/
void load_collmat_L (M, V, bpp, xm, ym, num, loc, start)
struct last_proc_collmat *M;
double *V;
int bpp, xm, ym, num, loc, start;
{
int ctr, n, m, bppm2, bppm1, locm2, locm1;

bppm1 = bpp - 1;
bppm2 = bpp - 2;
locm1 = loc - 1;
locm2 = loc - 2;
ctr = start - 1;

if (xm < ym)
    {
    if (num == 0)  /* xl_yg */
        {
	for (n = -1; n <= bppm2; n++)
	    {
	    unload_W_L_n0 (M, V, &ctr, n);
	    for (m = 1; m <= locm2; m++)
	        unload_int_L_n0 (M, V, &ctr, n, m);
	    unload_E_L_n0 (M, V, &ctr, n, locm1);
	    }
	unload_NW_L_n0 (M, V, &ctr, bppm1);
	for (m = 1; m <= locm2; m++)
	    unload_N_L_n0 (M, V, &ctr, bppm1, m);
	unload_NE_L_n0 (M, V, &ctr, bppm1, locm1);
        }

    else  /* num == 1, yl_xg */
        {
	for (m = -1; m <= bppm2; m++)
	    unload_S_L_n1 (M, V, &ctr, m);
	unload_SE_L_n1 (M, V, &ctr, bppm1);
	for (n = 1; n <= locm2; n++)
	    {
	    for (m = -1; m <= bppm2; m++)
	        unload_int_L_n1 (M, V, &ctr, n, m);
	    unload_E_L_n1 (M, V, &ctr, n, bppm1);
	    }
	for (m = -1; m <= bppm2; m++)
	    unload_N_L_n1 (M, V, &ctr, locm1, m);
	unload_NE_L_n1 (M, V, &ctr, locm1, bppm1);
        }
  }

else  /* xm >= ym */
    {
    if (num == 0)  /* xl_yg */
        {
	for (n = -1; n <= bppm2; n++)
	    unload_W_L_n0 (M, V, &ctr, n);
	unload_NW_L_n0 (M, V, &ctr, bppm1);
	for (m = 1; m <= locm2; m++)
	    {
	    for (n = -1; n <= bppm2; n++)
	        unload_int_L_n0 (M, V, &ctr, n, m);
	    unload_N_L_n0 (M, V, &ctr, bppm1, m);
	    }
	for (n = -1; n <= bppm2; n++)
	    unload_E_L_n0 (M, V, &ctr, n, locm1);
	unload_NE_L_n0 (M, V, &ctr, bppm1, locm1);
        }

    else  /* num == 1, yl_xg */
        {
	for (m = -1; m <= bppm2; m++)
	    {
	    unload_S_L_n1 (M, V, &ctr, m);
	    for (n = 1; n <= locm2; n++)
	        unload_int_L_n1 (M, V, &ctr, n, m);
	    unload_N_L_n1 (M, V, &ctr, locm1, m);
	    }
	unload_SE_L_n1 (M, V, &ctr, bppm1);
	for (n = 1; n <= locm2; n++)
	    unload_E_L_n1 (M, V, &ctr, n, bppm1);
	unload_NE_L_n1 (M, V, &ctr, locm1, bppm1);
        }
    }
}  /* end of function load_collmat_L */