454-calculix Report

**454-calculix**
Exec Time (%)	Codelet Name	Error (%)
66.4	CalculiX_main_181	42.7
3.6	mastruct_mastruct_615	0.9
9.2	isortii_isortii__199	18.7
3.8	SPOOLES_FrontMtx_src_FrontMtx_update_Fr...	85.7
3.8	SPOOLES_Chv_src_Chv_update_Chv_updateS_...	6.0
8.4	mafillsm_mafillsm__146	1.2
2.1	e_c3d_e_c3d__247	14.8
1.9	SPOOLES_Utilities_src_Utilities_newsort...	6.6
9.2	isortii_isortii__199	18.7
3.2	mastruct_mastruct_705	0.1
9.2	isortii_isortii__199	18.7
1.3	SPOOLES_Chv_src_Chv_factor_symmetric1x1...	38.9
14.8	mastruct_mastruct_331	1.0
2.5	SPOOLES_Utilities_src_Utilities_sortAnd...	1.3
3.5	results_results__1301	97.2
2.6	SPOOLES_InpMtx_src_InpMtx_fullAdj_InpMt...	5.9
1.6	results_results__1256	94.7
1.8	spooles_spooles_86	7.1
1.7	CalculiX_main_145	99.9

Exec Time : 1.306602e+08 cycles

Matching : 40.0%

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivoSPOOLES_Chv_src_Chv_factor_symmetric1x1_600

Real time: 1.753996e+06 -- Predicted time: 2.873963e+06

Callcount: 720

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
198	### 1	291.171626984127	5.787000e+03	4.032000e+03	30.3
362	### 2	505.937269372694	2.350000e+03	1.084000e+03	53.8

                /*     CalculiX - A 3-dimensional finite element program                 */
/*              Copyright (C) 1998 Guido Dhondt                          */

/*     This program is free software; you can redistribute it and/or     */
/*     modify it under the terms of the GNU General Public License as    */
/*     published by the Free Software Foundation(version 2);    */
/*                                                                       */

/*     This program is distributed in the hope that it will be useful,   */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
/*     GNU General Public License for more details.                      */

/*     You should have received a copy of the GNU General Public License */
/*     along with this program; if not, write to the Free Software       */
/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "CalculiX.h"

#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))

void mastruct(int *nk, int *kon, int *ipkon, char *lakon, int *ne,
	      int *nodeboun, int *ndirboun, int *nboun, int *ipompc,
	      int *nodempc, int *nmpc, int *nactdof, int *icol,
	      int *jq, int **mast1p, int **irowp, int *isolver, int *neq,
	      int *nnn, int *ikmpc, int *ilmpc, int *ikcol,
	      int *ipointer, int *nsky, int *nzs, int *nmethod,
              int *ithermal){

  int i,j,k,l,jj,ll,id,index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,
    ist1,ist2,node1,node2,isubtract,nmast,ifree,istart,istartold,itot,
    index1,index2,m,node,nzs_,ist,kflag,indexe,nope,isize,*mast1=NULL,
    *irow=NULL,nsky_exp,nsky_inc;

  /* the indices in the comments follow FORTRAN convention, i.e. the
     fields start with 1 */

  mast1=*mast1p;
  irow=*irowp;

  kflag=2;
  nzs_=*nzs;

  /* initialisation of nactmpc */

  for(i=0;i<4**nk;++i){nactdof[i]=0;}

  /* determining the mechanical active degrees of freedom due to elements */

  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      for(k=1;k<4;++k){
		  nactdof[4*node+k]=1;
	      }
	  }
      }
  }

  /* determining the thermal active degrees of freedom due to elements */

  if(*ithermal>1){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      nactdof[4*node]=1;
	  }
      }
  }

  /* determining the active degrees of freedom due to mpc's */

  for(i=0;i<*nmpc;++i){
      index=ipompc[i]-1;
      do{
	  nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=1;
	  index=nodempc[3*index+2];
	  if(index==0) break;
	  index--;
      }while(1);
  }
	   
  /* subtracting the SPC and MPC nodes */

  for(i=0;i<*nboun;++i){
      if(ndirboun[i]>3){continue;}
      nactdof[4*(nodeboun[i]-1)+ndirboun[i]]=0;
  }

  for(i=0;i<*nmpc;++i){
    index=ipompc[i]-1;
    nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=0;
  }
  
  /* numbering the active degrees of freedom */
  
  *neq=0;
  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*nk;++i){
	  for(j=1;j<4;++j){
	      if(nactdof[4*nnn[i]+j-4]!=0){
		  ++*neq;
		  nactdof[4*nnn[i]+j-4]=*neq;
	      }
	  }
      }
  }
  if(*ithermal>1){
      for(i=0;i<*nk;++i){
	  if(nactdof[4*nnn[i]-4]!=0){
	      ++*neq;
	      nactdof[4*nnn[i]-4]=*neq;
	  }
      }
  }
  
  ifree=0;
  
  /* determining the position of the skyline only */
  
  if(*isolver==1){
    
    /* ipointer(i) = first nonzero row number in column i */
    
    for(i=0;i<4**nk;++i){ipointer[i]=i+1;}
    
    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    ipointer[max(jdof1,jdof2)-1]=
	      min(ipointer[max(jdof1,jdof2)-1],min(jdof1,jdof2));
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    ipointer[max(idof1,idof2)-1]=
		      min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }
    
    /* defining icol and jq
       
       jq(i) = first nonzero row element in column i
       icol(i) = length of skyline up to and including column i */
    
    icol[0]=0;
    for(i=1;i<*neq;++i){
      jq[i]=ipointer[i];
      icol[i]=icol[i-1]+i-jq[i]+1;
    }
    *nsky=icol[*neq-1];

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("total length of the skyline\n");
    printf("%d\n\n",icol[*neq-1]);
    
  }
  else{

    /* determining the position of each nonzero matrix element

       mast1(ipointer(i)) = first nonzero row in column i
       irow(ipointer(i))  points to further nonzero elements in 
                             column i */
      
    for(i=0;i<4**nk;++i){ipointer[i]=0;}

    /* mechanical entries */
    
    if((*ithermal<2)||(*ithermal==3)){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }

    /* thermal entries*/

    if(*ithermal>1){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<nope;++jj){
	
	j=jj;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)];
	
	for(ll=jj;ll<nope;++ll){
	  
	  l=ll;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1-7;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2-7;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1-7;
	    idof2=7*node2-7;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }
    
    /* ordering the nonzero nodes in the SUPERdiagonal columns
       mast1 contains the row numbers column per column,
       irow the column numbers */

    for(i=0;i<*neq;++i){
      itot=0;
      if(ipointer[i]==0){
	printf("error in mastruct: zero column\n");
	FORTRAN(stop,());
      }
      istart=ipointer[i];
      while(1){
	++itot;
	ikcol[itot-1]=mast1[istart-1];
	istart=irow[istart-1];
	if(istart==0) break;
      }
      FORTRAN(isortii,(ikcol,icol,&itot,&kflag));
      istart=ipointer[i];
      for(j=0;j<itot-1;++j){
	mast1[istart-1]=ikcol[j];
	istartold=istart;
	istart=irow[istart-1];
	irow[istartold-1]=i+1;
      }
      mast1[istart-1]=ikcol[itot-1];
      irow[istart-1]=i+1;
    }

    /* definiing icol and jq */

    *nsky=0;
    nsky_exp=0;
    for(i=1;i<*neq;++i){
      nsky_inc=i+1-mast1[ipointer[i]-1];
      if(2147483647-*nsky<nsky_inc){
	++nsky_exp;
	nsky_inc=nsky_inc-2147483647;
      }
      *nsky=*nsky+nsky_inc;
    }

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("number of nonzero matrix elements\n");
    printf("%d\n",ifree);
    printf("total length of the skyline\n");
    printf("%d*2147483647+%d\n",nsky_exp,*nsky);
    printf("percentage of nonzero skyline elements\n");
    printf("%f\n\n",((double)ifree)/
	   ((double)(*nsky+*neq)+nsky_exp*(double)(2147483647))*100.);

    /* new meaning of icol,j1,mast1,irow:

       - irow is going to contain the row numbers of the SUBdiagonal
         nonzero's, column per column
       - mast1 contains the column numbers
       - icol(i)=# SUBdiagonal nonzero's in column i
       - jq(i)= location in field irow of the first SUBdiagonal
         nonzero in column i

	 */

    nmast=ifree;

    /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */

    FORTRAN(isortii,(mast1,irow,&nmast,&kflag));
    
    /* filtering out the diagonal elements and generating icol and jq */

    isubtract=0;
    for(i=1;i<*neq;++i){icol[i]=0;}
    k=0;
    for(i=0;i<nmast;++i){
      if(mast1[i]==irow[i]){++isubtract;}
      else{
	mast1[i-isubtract]=mast1[i];
	irow[i-isubtract]=irow[i];
	if(k!=mast1[i]){
	  for(l=k;l<mast1[i];++l){jq[l]=i+1-isubtract;}
	  k=mast1[i];
	}
	++icol[k-1];
      }
    }
    nmast=nmast-isubtract;
    for(l=k;l<*neq+1;++l){jq[l]=nmast+1;}

    for(i=0;i<*neq;++i){
      if(jq[i+1]-jq[i]>0){
	isize=jq[i+1]-jq[i];
	FORTRAN(isortii,(&irow[jq[i]-1],&mast1[jq[i]-1],&isize,&kflag));
      }
    }

    *nzs=jq[*neq-1]-1;

  }

  *mast1p=mast1;
  *irowp=irow;

  return;

}
	  
/*

What follows is the original FORTRAN code. The C Code is a one-to-one
manual translation of the FORTRAN code. However, the FORTRAN code might
be easier to understand.

!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine mastruct(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
     &  nboun,ipompc,
     &  nodempc,nmpc,nactdof,icol,jq,mast1,irow,isolver,neq,nnn,
     &  ikmpc,ilmpc,ikcol,ipointer,nsky,nzs,nmethod)
!
      implicit none
!
      character*6 lakon(*)
!
      integer kon(*),nodeboun(*),ndirboun(*),nodempc(3,*),ipompc(*),
     &  nactdof(3,*),icol(*),jq(*),ipointer(*),nnn(*),ikmpc(*),ilmpc(*),
     &  ikcol(*),mast1(*),irow(*),ipkon(*)
!
      integer nk,ne,nboun,nmpc,isolver,neq,nsky,nzs,i,j,k,l,jj,ll,id,
     &  index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,ist1,ist2,node1,
     &  node2,isubtract,nmast,ifree,istart,istartold,itot,index1,index2,
     &  m,node,nzs_,ist,kflag,nmethod,indexe,nope,nsky_exp,nsky_inc
!
      kflag=2
      nzs_=nzs
!
!     initialisation of nactmpc
!
      do i=1,nk
        do j=1,3
          nactdof(j,i)=0
        enddo
      enddo
!
!     determining the active degrees of freedom due to elements
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
            nope=20
         elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
            nope=8
         else
            nope=10
         endif
!
         do j=1,nope
            node=kon(indexe+j)
            do k=1,3
               nactdof(k,node)=1
            enddo
         enddo
      enddo
!
!     determining the active degrees of freedom due to mpc's
!
      do i=1,nmpc
         index=ipompc(i)
         do
            nactdof(nodempc(2,index),nodempc(1,index))=1
            index=nodempc(3,index)
	    if(index.eq.0) exit
	 enddo
      enddo
!      
!     subtracting the SPC and MPC nodes
!
      do i=1,nboun
        nactdof(ndirboun(i),nodeboun(i))=0
      enddo
!
      do i=1,nmpc
         index=ipompc(i)
         nactdof(nodempc(2,index),nodempc(1,index))=0
      enddo
!
!     numbering the active degrees of freedom
!
      neq=0
      do i=1,nk
        do j=1,3
          if(nactdof(j,nnn(i)).ne.0) then
            neq=neq+1
            nactdof(j,nnn(i))=neq
          endif
        enddo
      enddo
!
      ifree=0
!
!     determining the position of the skyline only
!
      if(isolver.eq.1) then
!
!       ipointer(i) = first nonzero row number in column i
!
        do i=1,3*nk
          ipointer(i)=i
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 ipointer(max(jdof1,jdof2))=
     &             min(ipointer(max(jdof1,jdof2)),min(jdof1,jdof2))
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       defining icol and jq
!
!       jq(i) = first nonzero row element in column i
!       icol(i) = length of skyline up to and including column i
!
        icol(1)=0
        do i=2,neq
          jq(i)=ipointer(i)
          icol(i)=icol(i-1)+i-jq(i)
        enddo
        nsky=icol(neq)
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'total length of the skyline'
        write(*,*) icol(neq)
        write(*,*)
!
      else
!
!     determining the position of each nonzero matrix element
!
!       mast1(ipointer(i)) = first nonzero row in column i
!       irow(ipointer(i)) points to further nonzero elements in
!         column i
!
        do i=1,3*nk
          ipointer(i)=0
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 call insert(ipointer,mast1,irow,
     &                       jdof1,jdof2,ifree,nzs_)
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                          call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       ordering the nonzero nodes in the SUPERdiagonal columns
!       mast1 contains the row numbers column per column,
!       irow the column numbers
!
        do i=1,neq
          itot=0
          if(ipointer(i).eq.0) then
            write(*,*) 'error in mastruct: zero column'
            stop
          endif
          istart=ipointer(i)
          do
            itot=itot+1
            ikcol(itot)=mast1(istart)
            istart=irow(istart)
            if(istart.eq.0) exit
          enddo
          call isortii(ikcol,icol,itot,kflag)
          istart=ipointer(i)
          do j=1,itot-1
            mast1(istart)=ikcol(j)
            istartold=istart
            istart=irow(istart)
            irow(istartold)=i
          enddo
          mast1(istart)=ikcol(itot)
          irow(istart)=i
        enddo
!
!       defining icol and jq
!
        nsky=0
	nsky_exp=0
        do i=2,neq
          nsky_inc=i-mast1(ipointer(i))
          if(2147483647-nsky.lt.nsky_inc) then
	     nsky_exp=nsky_exp+1
	     nsky_inc=nsky_inc-2147483647
	  endif
	  nsky=nsky+nsky_inc
        enddo
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'number of nonzero matrix elements'
        write(*,*) ifree
        write(*,*) 'total length of the skyline'
        write(*,*) nsky_exp,'*2147483647+',nsky
        write(*,*) 'percentage of nonzero skyline elements'
        write(*,*) real(ifree)/
     &     (real(nsky+neq)+nsky_exp*real(2147483647))*100.
        write(*,*)
!
!       new meaning of icol,j1,mast1,irow:
!       - irow is going to contain the row numbers of the SUBdiagonal
!         nonzero's, column per column
!       - mast1 contains the column numbers
!       - icol(i)=# SUBdiagonal nonzero's in column i
!       - jq(i)= location in field irow of the first SUBdiagonal
!         nonzero in column i
!
        nmast=ifree
!
!       switching from a SUPERdiagonal inventary to a SUBdiagonal one
!
        call isortii(mast1,irow,nmast,kflag)
!
!       filtering out the diagonal elements and generating icol and jq
!
        isubtract=0
        do i=1,neq
           icol(i)=0
        enddo
        k=0
        do i=1,nmast
           if(mast1(i).eq.irow(i)) then
              isubtract=isubtract+1
           else
              mast1(i-isubtract)=mast1(i)
              irow(i-isubtract)=irow(i)
              if(k.ne.mast1(i)) then
                 do l=k+1,mast1(i)
                    jq(l)=i-isubtract
                 enddo
                 k=mast1(i)
              endif
              icol(k)=icol(k)+1
           endif
        enddo
        nmast=nmast-isubtract
        do l=k+1,neq+1
           jq(l)=nmast+1
        enddo
!
        do i=1,neq
           if(jq(i+1)-jq(i).gt.0) then
              call isortii(irow(jq(i)),mast1(jq(i)),jq(i+1)-jq(i),
     &          kflag)

           endif
        enddo
!
        nzs=jq(neq)-1
!
      endif
!
      return
      end

      */

invivomastruct_mastruct_615

Real time: 4.825128e+06 -- Predicted time: 4.780256e+06

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	4.780256e+06	4.786612e+06	0.1

                /*     CalculiX - A 3-dimensional finite element program                 */
/*              Copyright (C) 1998 Guido Dhondt                          */

/*     This program is free software; you can redistribute it and/or     */
/*     modify it under the terms of the GNU General Public License as    */
/*     published by the Free Software Foundation(version 2);    */
/*                                                                       */

/*     This program is distributed in the hope that it will be useful,   */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
/*     GNU General Public License for more details.                      */

/*     You should have received a copy of the GNU General Public License */
/*     along with this program; if not, write to the Free Software       */
/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "CalculiX.h"

#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))

void mastruct(int *nk, int *kon, int *ipkon, char *lakon, int *ne,
	      int *nodeboun, int *ndirboun, int *nboun, int *ipompc,
	      int *nodempc, int *nmpc, int *nactdof, int *icol,
	      int *jq, int **mast1p, int **irowp, int *isolver, int *neq,
	      int *nnn, int *ikmpc, int *ilmpc, int *ikcol,
	      int *ipointer, int *nsky, int *nzs, int *nmethod,
              int *ithermal){

  int i,j,k,l,jj,ll,id,index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,
    ist1,ist2,node1,node2,isubtract,nmast,ifree,istart,istartold,itot,
    index1,index2,m,node,nzs_,ist,kflag,indexe,nope,isize,*mast1=NULL,
    *irow=NULL,nsky_exp,nsky_inc;

  /* the indices in the comments follow FORTRAN convention, i.e. the
     fields start with 1 */

  mast1=*mast1p;
  irow=*irowp;

  kflag=2;
  nzs_=*nzs;

  /* initialisation of nactmpc */

  for(i=0;i<4**nk;++i){nactdof[i]=0;}

  /* determining the mechanical active degrees of freedom due to elements */

  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      for(k=1;k<4;++k){
		  nactdof[4*node+k]=1;
	      }
	  }
      }
  }

  /* determining the thermal active degrees of freedom due to elements */

  if(*ithermal>1){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      nactdof[4*node]=1;
	  }
      }
  }

  /* determining the active degrees of freedom due to mpc's */

  for(i=0;i<*nmpc;++i){
      index=ipompc[i]-1;
      do{
	  nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=1;
	  index=nodempc[3*index+2];
	  if(index==0) break;
	  index--;
      }while(1);
  }
	   
  /* subtracting the SPC and MPC nodes */

  for(i=0;i<*nboun;++i){
      if(ndirboun[i]>3){continue;}
      nactdof[4*(nodeboun[i]-1)+ndirboun[i]]=0;
  }

  for(i=0;i<*nmpc;++i){
    index=ipompc[i]-1;
    nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=0;
  }
  
  /* numbering the active degrees of freedom */
  
  *neq=0;
  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*nk;++i){
	  for(j=1;j<4;++j){
	      if(nactdof[4*nnn[i]+j-4]!=0){
		  ++*neq;
		  nactdof[4*nnn[i]+j-4]=*neq;
	      }
	  }
      }
  }
  if(*ithermal>1){
      for(i=0;i<*nk;++i){
	  if(nactdof[4*nnn[i]-4]!=0){
	      ++*neq;
	      nactdof[4*nnn[i]-4]=*neq;
	  }
      }
  }
  
  ifree=0;
  
  /* determining the position of the skyline only */
  
  if(*isolver==1){
    
    /* ipointer(i) = first nonzero row number in column i */
    
    for(i=0;i<4**nk;++i){ipointer[i]=i+1;}
    
    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    ipointer[max(jdof1,jdof2)-1]=
	      min(ipointer[max(jdof1,jdof2)-1],min(jdof1,jdof2));
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    ipointer[max(idof1,idof2)-1]=
		      min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }
    
    /* defining icol and jq
       
       jq(i) = first nonzero row element in column i
       icol(i) = length of skyline up to and including column i */
    
    icol[0]=0;
    for(i=1;i<*neq;++i){
      jq[i]=ipointer[i];
      icol[i]=icol[i-1]+i-jq[i]+1;
    }
    *nsky=icol[*neq-1];

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("total length of the skyline\n");
    printf("%d\n\n",icol[*neq-1]);
    
  }
  else{

    /* determining the position of each nonzero matrix element

       mast1(ipointer(i)) = first nonzero row in column i
       irow(ipointer(i))  points to further nonzero elements in 
                             column i */
      
    for(i=0;i<4**nk;++i){ipointer[i]=0;}

    /* mechanical entries */
    
    if((*ithermal<2)||(*ithermal==3)){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }

    /* thermal entries*/

    if(*ithermal>1){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<nope;++jj){
	
	j=jj;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)];
	
	for(ll=jj;ll<nope;++ll){
	  
	  l=ll;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1-7;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2-7;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1-7;
	    idof2=7*node2-7;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }
    
    /* ordering the nonzero nodes in the SUPERdiagonal columns
       mast1 contains the row numbers column per column,
       irow the column numbers */

    for(i=0;i<*neq;++i){
      itot=0;
      if(ipointer[i]==0){
	printf("error in mastruct: zero column\n");
	FORTRAN(stop,());
      }
      istart=ipointer[i];
      while(1){
	++itot;
	ikcol[itot-1]=mast1[istart-1];
	istart=irow[istart-1];
	if(istart==0) break;
      }
      FORTRAN(isortii,(ikcol,icol,&itot,&kflag));
      istart=ipointer[i];
      for(j=0;j<itot-1;++j){
	mast1[istart-1]=ikcol[j];
	istartold=istart;
	istart=irow[istart-1];
	irow[istartold-1]=i+1;
      }
      mast1[istart-1]=ikcol[itot-1];
      irow[istart-1]=i+1;
    }

    /* definiing icol and jq */

    *nsky=0;
    nsky_exp=0;
    for(i=1;i<*neq;++i){
      nsky_inc=i+1-mast1[ipointer[i]-1];
      if(2147483647-*nsky<nsky_inc){
	++nsky_exp;
	nsky_inc=nsky_inc-2147483647;
      }
      *nsky=*nsky+nsky_inc;
    }

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("number of nonzero matrix elements\n");
    printf("%d\n",ifree);
    printf("total length of the skyline\n");
    printf("%d*2147483647+%d\n",nsky_exp,*nsky);
    printf("percentage of nonzero skyline elements\n");
    printf("%f\n\n",((double)ifree)/
	   ((double)(*nsky+*neq)+nsky_exp*(double)(2147483647))*100.);

    /* new meaning of icol,j1,mast1,irow:

       - irow is going to contain the row numbers of the SUBdiagonal
         nonzero's, column per column
       - mast1 contains the column numbers
       - icol(i)=# SUBdiagonal nonzero's in column i
       - jq(i)= location in field irow of the first SUBdiagonal
         nonzero in column i

	 */

    nmast=ifree;

    /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */

    FORTRAN(isortii,(mast1,irow,&nmast,&kflag));
    
    /* filtering out the diagonal elements and generating icol and jq */

    isubtract=0;
    for(i=1;i<*neq;++i){icol[i]=0;}
    k=0;
    for(i=0;i<nmast;++i){
      if(mast1[i]==irow[i]){++isubtract;}
      else{
	mast1[i-isubtract]=mast1[i];
	irow[i-isubtract]=irow[i];
	if(k!=mast1[i]){
	  for(l=k;l<mast1[i];++l){jq[l]=i+1-isubtract;}
	  k=mast1[i];
	}
	++icol[k-1];
      }
    }
    nmast=nmast-isubtract;
    for(l=k;l<*neq+1;++l){jq[l]=nmast+1;}

    for(i=0;i<*neq;++i){
      if(jq[i+1]-jq[i]>0){
	isize=jq[i+1]-jq[i];
	FORTRAN(isortii,(&irow[jq[i]-1],&mast1[jq[i]-1],&isize,&kflag));
      }
    }

    *nzs=jq[*neq-1]-1;

  }

  *mast1p=mast1;
  *irowp=irow;

  return;

}
	  
/*

What follows is the original FORTRAN code. The C Code is a one-to-one
manual translation of the FORTRAN code. However, the FORTRAN code might
be easier to understand.

!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine mastruct(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
     &  nboun,ipompc,
     &  nodempc,nmpc,nactdof,icol,jq,mast1,irow,isolver,neq,nnn,
     &  ikmpc,ilmpc,ikcol,ipointer,nsky,nzs,nmethod)
!
      implicit none
!
      character*6 lakon(*)
!
      integer kon(*),nodeboun(*),ndirboun(*),nodempc(3,*),ipompc(*),
     &  nactdof(3,*),icol(*),jq(*),ipointer(*),nnn(*),ikmpc(*),ilmpc(*),
     &  ikcol(*),mast1(*),irow(*),ipkon(*)
!
      integer nk,ne,nboun,nmpc,isolver,neq,nsky,nzs,i,j,k,l,jj,ll,id,
     &  index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,ist1,ist2,node1,
     &  node2,isubtract,nmast,ifree,istart,istartold,itot,index1,index2,
     &  m,node,nzs_,ist,kflag,nmethod,indexe,nope,nsky_exp,nsky_inc
!
      kflag=2
      nzs_=nzs
!
!     initialisation of nactmpc
!
      do i=1,nk
        do j=1,3
          nactdof(j,i)=0
        enddo
      enddo
!
!     determining the active degrees of freedom due to elements
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
            nope=20
         elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
            nope=8
         else
            nope=10
         endif
!
         do j=1,nope
            node=kon(indexe+j)
            do k=1,3
               nactdof(k,node)=1
            enddo
         enddo
      enddo
!
!     determining the active degrees of freedom due to mpc's
!
      do i=1,nmpc
         index=ipompc(i)
         do
            nactdof(nodempc(2,index),nodempc(1,index))=1
            index=nodempc(3,index)
	    if(index.eq.0) exit
	 enddo
      enddo
!      
!     subtracting the SPC and MPC nodes
!
      do i=1,nboun
        nactdof(ndirboun(i),nodeboun(i))=0
      enddo
!
      do i=1,nmpc
         index=ipompc(i)
         nactdof(nodempc(2,index),nodempc(1,index))=0
      enddo
!
!     numbering the active degrees of freedom
!
      neq=0
      do i=1,nk
        do j=1,3
          if(nactdof(j,nnn(i)).ne.0) then
            neq=neq+1
            nactdof(j,nnn(i))=neq
          endif
        enddo
      enddo
!
      ifree=0
!
!     determining the position of the skyline only
!
      if(isolver.eq.1) then
!
!       ipointer(i) = first nonzero row number in column i
!
        do i=1,3*nk
          ipointer(i)=i
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 ipointer(max(jdof1,jdof2))=
     &             min(ipointer(max(jdof1,jdof2)),min(jdof1,jdof2))
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       defining icol and jq
!
!       jq(i) = first nonzero row element in column i
!       icol(i) = length of skyline up to and including column i
!
        icol(1)=0
        do i=2,neq
          jq(i)=ipointer(i)
          icol(i)=icol(i-1)+i-jq(i)
        enddo
        nsky=icol(neq)
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'total length of the skyline'
        write(*,*) icol(neq)
        write(*,*)
!
      else
!
!     determining the position of each nonzero matrix element
!
!       mast1(ipointer(i)) = first nonzero row in column i
!       irow(ipointer(i)) points to further nonzero elements in
!         column i
!
        do i=1,3*nk
          ipointer(i)=0
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 call insert(ipointer,mast1,irow,
     &                       jdof1,jdof2,ifree,nzs_)
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                          call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       ordering the nonzero nodes in the SUPERdiagonal columns
!       mast1 contains the row numbers column per column,
!       irow the column numbers
!
        do i=1,neq
          itot=0
          if(ipointer(i).eq.0) then
            write(*,*) 'error in mastruct: zero column'
            stop
          endif
          istart=ipointer(i)
          do
            itot=itot+1
            ikcol(itot)=mast1(istart)
            istart=irow(istart)
            if(istart.eq.0) exit
          enddo
          call isortii(ikcol,icol,itot,kflag)
          istart=ipointer(i)
          do j=1,itot-1
            mast1(istart)=ikcol(j)
            istartold=istart
            istart=irow(istart)
            irow(istartold)=i
          enddo
          mast1(istart)=ikcol(itot)
          irow(istart)=i
        enddo
!
!       defining icol and jq
!
        nsky=0
	nsky_exp=0
        do i=2,neq
          nsky_inc=i-mast1(ipointer(i))
          if(2147483647-nsky.lt.nsky_inc) then
	     nsky_exp=nsky_exp+1
	     nsky_inc=nsky_inc-2147483647
	  endif
	  nsky=nsky+nsky_inc
        enddo
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'number of nonzero matrix elements'
        write(*,*) ifree
        write(*,*) 'total length of the skyline'
        write(*,*) nsky_exp,'*2147483647+',nsky
        write(*,*) 'percentage of nonzero skyline elements'
        write(*,*) real(ifree)/
     &     (real(nsky+neq)+nsky_exp*real(2147483647))*100.
        write(*,*)
!
!       new meaning of icol,j1,mast1,irow:
!       - irow is going to contain the row numbers of the SUBdiagonal
!         nonzero's, column per column
!       - mast1 contains the column numbers
!       - icol(i)=# SUBdiagonal nonzero's in column i
!       - jq(i)= location in field irow of the first SUBdiagonal
!         nonzero in column i
!
        nmast=ifree
!
!       switching from a SUPERdiagonal inventary to a SUBdiagonal one
!
        call isortii(mast1,irow,nmast,kflag)
!
!       filtering out the diagonal elements and generating icol and jq
!
        isubtract=0
        do i=1,neq
           icol(i)=0
        enddo
        k=0
        do i=1,nmast
           if(mast1(i).eq.irow(i)) then
              isubtract=isubtract+1
           else
              mast1(i-isubtract)=mast1(i)
              irow(i-isubtract)=irow(i)
              if(k.ne.mast1(i)) then
                 do l=k+1,mast1(i)
                    jq(l)=i-isubtract
                 enddo
                 k=mast1(i)
              endif
              icol(k)=icol(k)+1
           endif
        enddo
        nmast=nmast-isubtract
        do l=k+1,neq+1
           jq(l)=nmast+1
        enddo
!
        do i=1,neq
           if(jq(i+1)-jq(i).gt.0) then
              call isortii(irow(jq(i)),mast1(jq(i)),jq(i+1)-jq(i),
     &          kflag)

           endif
        enddo
!
        nzs=jq(neq)-1
!
      endif
!
      return
      end

      */

invivomastruct_mastruct_331

Real time: 1.939523e+07 -- Predicted time: 1.919306e+07

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	1.919306e+07	1.953532e+07	1.7

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivoSPOOLES_Utilities_src_Utilities_DV_DVdot33_1241

Real time: 3.388256e+06 -- Predicted time: 7.024252e+06

Callcount: 7920

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1054	### 1	1004.93464052288	1.065000e+03	6.120000e+02	42.5
6	### 2	1025.04587155963	9.890000e+02	4.360000e+02	55.9
14	### 3	446.889830508475	1.072000e+03	4.720000e+02	55.9
7	### 4	693.803571428571	8.630000e+02	4.480000e+02	48.0
644	### 5	812.237037037037	8.410000e+02	5.400000e+02	35.7
2944	### 6	410.880681818182	1.030000e+03	7.040000e+02	31.6
1410	### 7	1474.89655172414	6.930000e+02	3.480000e+02	49.7
1308	### 8	1410.60784313725	7.960000e+02	2.040000e+02	74.3
5547	### 9	670.869791666667	9.110000e+02	7.680000e+02	15.6

                /*     CalculiX - A 3-dimensional finite element program                 */
/*              Copyright (C) 1998 Guido Dhondt                          */

/*     This program is free software; you can redistribute it and/or     */
/*     modify it under the terms of the GNU General Public License as    */
/*     published by the Free Software Foundation(version 2);    */
/*                                                                       */

/*     This program is distributed in the hope that it will be useful,   */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
/*     GNU General Public License for more details.                      */

/*     You should have received a copy of the GNU General Public License */
/*     along with this program; if not, write to the Free Software       */
/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "CalculiX.h"

#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))

void mastruct(int *nk, int *kon, int *ipkon, char *lakon, int *ne,
	      int *nodeboun, int *ndirboun, int *nboun, int *ipompc,
	      int *nodempc, int *nmpc, int *nactdof, int *icol,
	      int *jq, int **mast1p, int **irowp, int *isolver, int *neq,
	      int *nnn, int *ikmpc, int *ilmpc, int *ikcol,
	      int *ipointer, int *nsky, int *nzs, int *nmethod,
              int *ithermal){

  int i,j,k,l,jj,ll,id,index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,
    ist1,ist2,node1,node2,isubtract,nmast,ifree,istart,istartold,itot,
    index1,index2,m,node,nzs_,ist,kflag,indexe,nope,isize,*mast1=NULL,
    *irow=NULL,nsky_exp,nsky_inc;

  /* the indices in the comments follow FORTRAN convention, i.e. the
     fields start with 1 */

  mast1=*mast1p;
  irow=*irowp;

  kflag=2;
  nzs_=*nzs;

  /* initialisation of nactmpc */

  for(i=0;i<4**nk;++i){nactdof[i]=0;}

  /* determining the mechanical active degrees of freedom due to elements */

  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      for(k=1;k<4;++k){
		  nactdof[4*node+k]=1;
	      }
	  }
      }
  }

  /* determining the thermal active degrees of freedom due to elements */

  if(*ithermal>1){
      for(i=0;i<*ne;++i){
	  
	  if(ipkon[i]<0) continue;
	  indexe=ipkon[i];
	  if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
	  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
	  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
	  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
	  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
	  else {nope=6;}
	  
	  for(j=0;j<nope;++j){
	      node=kon[indexe+j]-1;
	      nactdof[4*node]=1;
	  }
      }
  }

  /* determining the active degrees of freedom due to mpc's */

  for(i=0;i<*nmpc;++i){
      index=ipompc[i]-1;
      do{
	  nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=1;
	  index=nodempc[3*index+2];
	  if(index==0) break;
	  index--;
      }while(1);
  }
	   
  /* subtracting the SPC and MPC nodes */

  for(i=0;i<*nboun;++i){
      if(ndirboun[i]>3){continue;}
      nactdof[4*(nodeboun[i]-1)+ndirboun[i]]=0;
  }

  for(i=0;i<*nmpc;++i){
    index=ipompc[i]-1;
    nactdof[4*nodempc[3*index]+nodempc[3*index+1]-4]=0;
  }
  
  /* numbering the active degrees of freedom */
  
  *neq=0;
  if((*ithermal<2)||(*ithermal==3)){
      for(i=0;i<*nk;++i){
	  for(j=1;j<4;++j){
	      if(nactdof[4*nnn[i]+j-4]!=0){
		  ++*neq;
		  nactdof[4*nnn[i]+j-4]=*neq;
	      }
	  }
      }
  }
  if(*ithermal>1){
      for(i=0;i<*nk;++i){
	  if(nactdof[4*nnn[i]-4]!=0){
	      ++*neq;
	      nactdof[4*nnn[i]-4]=*neq;
	  }
      }
  }
  
  ifree=0;
  
  /* determining the position of the skyline only */
  
  if(*isolver==1){
    
    /* ipointer(i) = first nonzero row number in column i */
    
    for(i=0;i<4**nk;++i){ipointer[i]=i+1;}
    
    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    ipointer[max(jdof1,jdof2)-1]=
	      min(ipointer[max(jdof1,jdof2)-1],min(jdof1,jdof2));
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    ipointer[max(idof1,idof2)-1]=
		      min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      ipointer[max(idof1,idof2)-1]=
			min(ipointer[max(idof1,idof2)-1],min(idof1,idof2));}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }
    
    /* defining icol and jq
       
       jq(i) = first nonzero row element in column i
       icol(i) = length of skyline up to and including column i */
    
    icol[0]=0;
    for(i=1;i<*neq;++i){
      jq[i]=ipointer[i];
      icol[i]=icol[i-1]+i-jq[i]+1;
    }
    *nsky=icol[*neq-1];

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("total length of the skyline\n");
    printf("%d\n\n",icol[*neq-1]);
    
  }
  else{

    /* determining the position of each nonzero matrix element

       mast1(ipointer(i)) = first nonzero row in column i
       irow(ipointer(i))  points to further nonzero elements in 
                             column i */
      
    for(i=0;i<4**nk;++i){ipointer[i]=0;}

    /* mechanical entries */
    
    if((*ithermal<2)||(*ithermal==3)){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<3*nope;++jj){
	
	j=jj/3;
	k=jj-3*j;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)+k+1];
	
	for(ll=jj;ll<3*nope;++ll){
	  
	  l=ll/3;
	  m=ll-3*l;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)+m+1];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1+k-6;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2+m-6;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1+k-6;
	    idof2=7*node2+m-6;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }

    /* thermal entries*/

    if(*ithermal>1){

    for(i=0;i<*ne;++i){
      
      if(ipkon[i]<0) continue;
      indexe=ipkon[i];
      if(strcmp1(&lakon[8*i+3],"2")==0)nope=20;
      else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
      else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
      else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
      else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
      else {nope=6;}
      
      for(jj=0;jj<nope;++jj){
	
	j=jj;
	
	node1=kon[indexe+j];
	jdof1=nactdof[4*(node1-1)];
	
	for(ll=jj;ll<nope;++ll){
	  
	  l=ll;
	  
	  node2=kon[indexe+l];
	  jdof2=nactdof[4*(node2-1)];
	  
	  /* check whether one of the DOF belongs to a SPC or MPC */
	  
	  if((jdof1!=0)&&(jdof2!=0)){
	    insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
	  }
	  else if((jdof1!=0)||(jdof2!=0)){
	    
	    /* idof1: genuine DOF
	       idof2: nominal DOF of the SPC/MPC */
	    
	    if(jdof1==0){
	      idof1=jdof2;
	      idof2=7*node1-7;}
	    else{
	      idof1=jdof1;
	      idof2=7*node2-7;}
	    
	    if(*nmpc>0){
	      
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
	      if((id>0)&&(ikmpc[id-1]==idof2)){
		
		/* regular DOF / MPC */
		
		id=ilmpc[id-1];
		ist=ipompc[id-1];
		index=nodempc[3*ist-1];
		if(index==0) continue;
		while(1){
		  idof2=nactdof[4*nodempc[3*index-3]+nodempc[3*index-2]-4];
		  if(idof2!=0){
		    insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
		  }
		  index=nodempc[3*index-1];
		  if(index==0) break;
		}
		continue;
	      }
	    }
	  }
	  
	  else{
	    idof1=7*node1-7;
	    idof2=7*node2-7;
	    mpc1=0;
	    mpc2=0;
	    if(*nmpc>0){
	      FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
	      if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
	      FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
	      if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
	    }
	    if((mpc1==1)&&(mpc2==1)){
	      id1=ilmpc[id1-1];
	      id2=ilmpc[id2-1];
	      if(id1==id2){
		
		/* MPC id1 / MPC id1 */
		
		ist=ipompc[id1-1];
		index1=nodempc[3*ist-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  index2=index1;
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	      
	      else{
		
		/* MPC id1 /MPC id2 */
		
		ist1=ipompc[id1-1];
		index1=nodempc[3*ist1-1];
		if(index1==0) continue;
		while(1){
		  idof1=nactdof[4*nodempc[3*index1-3]+nodempc[3*index1-2]-4];
		  ist2=ipompc[id2-1];
		  index2=nodempc[3*ist2-1];
		  if(index2==0){
		    index1=nodempc[3*index1-1];
		    if(index1==0){break;}
		    else{continue;}
		  }
		  while(1){
		    idof2=nactdof[4*nodempc[3*index2-3]+nodempc[3*index2-2]-4];
		    if((idof1!=0)&&(idof2!=0)){
		      insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
		    index2=nodempc[3*index2-1];
		    if(index2==0) break;
		  }
		  index1=nodempc[3*index1-1];
		  if(index1==0) break;
		}
	      }
	    }
	  }
	}
      }
    }

    }
    
    /* ordering the nonzero nodes in the SUPERdiagonal columns
       mast1 contains the row numbers column per column,
       irow the column numbers */

    for(i=0;i<*neq;++i){
      itot=0;
      if(ipointer[i]==0){
	printf("error in mastruct: zero column\n");
	FORTRAN(stop,());
      }
      istart=ipointer[i];
      while(1){
	++itot;
	ikcol[itot-1]=mast1[istart-1];
	istart=irow[istart-1];
	if(istart==0) break;
      }
      FORTRAN(isortii,(ikcol,icol,&itot,&kflag));
      istart=ipointer[i];
      for(j=0;j<itot-1;++j){
	mast1[istart-1]=ikcol[j];
	istartold=istart;
	istart=irow[istart-1];
	irow[istartold-1]=i+1;
      }
      mast1[istart-1]=ikcol[itot-1];
      irow[istart-1]=i+1;
    }

    /* definiing icol and jq */

    *nsky=0;
    nsky_exp=0;
    for(i=1;i<*neq;++i){
      nsky_inc=i+1-mast1[ipointer[i]-1];
      if(2147483647-*nsky<nsky_inc){
	++nsky_exp;
	nsky_inc=nsky_inc-2147483647;
      }
      *nsky=*nsky+nsky_inc;
    }

    if(*neq==0){
      printf("*ERROR: no degrees of freedom in the model\n");
      FORTRAN(stop,());
    }

    printf("number of equations\n");
    printf("%d\n",*neq);
    printf("number of nonzero matrix elements\n");
    printf("%d\n",ifree);
    printf("total length of the skyline\n");
    printf("%d*2147483647+%d\n",nsky_exp,*nsky);
    printf("percentage of nonzero skyline elements\n");
    printf("%f\n\n",((double)ifree)/
	   ((double)(*nsky+*neq)+nsky_exp*(double)(2147483647))*100.);

    /* new meaning of icol,j1,mast1,irow:

       - irow is going to contain the row numbers of the SUBdiagonal
         nonzero's, column per column
       - mast1 contains the column numbers
       - icol(i)=# SUBdiagonal nonzero's in column i
       - jq(i)= location in field irow of the first SUBdiagonal
         nonzero in column i

	 */

    nmast=ifree;

    /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */

    FORTRAN(isortii,(mast1,irow,&nmast,&kflag));
    
    /* filtering out the diagonal elements and generating icol and jq */

    isubtract=0;
    for(i=1;i<*neq;++i){icol[i]=0;}
    k=0;
    for(i=0;i<nmast;++i){
      if(mast1[i]==irow[i]){++isubtract;}
      else{
	mast1[i-isubtract]=mast1[i];
	irow[i-isubtract]=irow[i];
	if(k!=mast1[i]){
	  for(l=k;l<mast1[i];++l){jq[l]=i+1-isubtract;}
	  k=mast1[i];
	}
	++icol[k-1];
      }
    }
    nmast=nmast-isubtract;
    for(l=k;l<*neq+1;++l){jq[l]=nmast+1;}

    for(i=0;i<*neq;++i){
      if(jq[i+1]-jq[i]>0){
	isize=jq[i+1]-jq[i];
	FORTRAN(isortii,(&irow[jq[i]-1],&mast1[jq[i]-1],&isize,&kflag));
      }
    }

    *nzs=jq[*neq-1]-1;

  }

  *mast1p=mast1;
  *irowp=irow;

  return;

}
	  
/*

What follows is the original FORTRAN code. The C Code is a one-to-one
manual translation of the FORTRAN code. However, the FORTRAN code might
be easier to understand.

!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine mastruct(nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
     &  nboun,ipompc,
     &  nodempc,nmpc,nactdof,icol,jq,mast1,irow,isolver,neq,nnn,
     &  ikmpc,ilmpc,ikcol,ipointer,nsky,nzs,nmethod)
!
      implicit none
!
      character*6 lakon(*)
!
      integer kon(*),nodeboun(*),ndirboun(*),nodempc(3,*),ipompc(*),
     &  nactdof(3,*),icol(*),jq(*),ipointer(*),nnn(*),ikmpc(*),ilmpc(*),
     &  ikcol(*),mast1(*),irow(*),ipkon(*)
!
      integer nk,ne,nboun,nmpc,isolver,neq,nsky,nzs,i,j,k,l,jj,ll,id,
     &  index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,ist1,ist2,node1,
     &  node2,isubtract,nmast,ifree,istart,istartold,itot,index1,index2,
     &  m,node,nzs_,ist,kflag,nmethod,indexe,nope,nsky_exp,nsky_inc
!
      kflag=2
      nzs_=nzs
!
!     initialisation of nactmpc
!
      do i=1,nk
        do j=1,3
          nactdof(j,i)=0
        enddo
      enddo
!
!     determining the active degrees of freedom due to elements
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
            nope=20
         elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
            nope=8
         else
            nope=10
         endif
!
         do j=1,nope
            node=kon(indexe+j)
            do k=1,3
               nactdof(k,node)=1
            enddo
         enddo
      enddo
!
!     determining the active degrees of freedom due to mpc's
!
      do i=1,nmpc
         index=ipompc(i)
         do
            nactdof(nodempc(2,index),nodempc(1,index))=1
            index=nodempc(3,index)
	    if(index.eq.0) exit
	 enddo
      enddo
!      
!     subtracting the SPC and MPC nodes
!
      do i=1,nboun
        nactdof(ndirboun(i),nodeboun(i))=0
      enddo
!
      do i=1,nmpc
         index=ipompc(i)
         nactdof(nodempc(2,index),nodempc(1,index))=0
      enddo
!
!     numbering the active degrees of freedom
!
      neq=0
      do i=1,nk
        do j=1,3
          if(nactdof(j,nnn(i)).ne.0) then
            neq=neq+1
            nactdof(j,nnn(i))=neq
          endif
        enddo
      enddo
!
      ifree=0
!
!     determining the position of the skyline only
!
      if(isolver.eq.1) then
!
!       ipointer(i) = first nonzero row number in column i
!
        do i=1,3*nk
          ipointer(i)=i
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 ipointer(max(jdof1,jdof2))=
     &             min(ipointer(max(jdof1,jdof2)),min(jdof1,jdof2))
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         ipointer(max(idof1,idof2))=
     &                         min(ipointer(max(idof1,idof2)),
     &                             min(idof1,idof2))
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       defining icol and jq
!
!       jq(i) = first nonzero row element in column i
!       icol(i) = length of skyline up to and including column i
!
        icol(1)=0
        do i=2,neq
          jq(i)=ipointer(i)
          icol(i)=icol(i-1)+i-jq(i)
        enddo
        nsky=icol(neq)
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'total length of the skyline'
        write(*,*) icol(neq)
        write(*,*)
!
      else
!
!     determining the position of each nonzero matrix element
!
!       mast1(ipointer(i)) = first nonzero row in column i
!       irow(ipointer(i)) points to further nonzero elements in
!         column i
!
        do i=1,3*nk
          ipointer(i)=0
        enddo
!
        do i=1,ne
!
          if(ipkon(i).lt.0) cycle
          indexe=ipkon(i)
          if((lakon(i).eq.'C3D20R').or.(lakon(i).eq.'C3D20 ')) then
             nope=20
          elseif((lakon(i).eq.'C3D8R ').or.(lakon(i).eq.'C3D8  ')) then
             nope=8
          else
             nope=10
          endif
!
          do jj=1,3*nope
!
            j=(jj-1)/3+1
            k=jj-3*(j-1)
!
            node1=kon(indexe+j)
            jdof1=nactdof(k,node1)
!
            do ll=jj,3*nope
!
              l=(ll-1)/3+1
              m=ll-3*(l-1)
!
              node2=kon(indexe+l)
              jdof2=nactdof(m,node2)
!
!             check whether one of the DOF belongs to a SPC or MPC
!
              if((jdof1.ne.0).and.(jdof2.ne.0)) then
                 call insert(ipointer,mast1,irow,
     &                       jdof1,jdof2,ifree,nzs_)
              elseif((jdof1.ne.0).or.(jdof2.ne.0)) then
!
!              idof1: genuine DOF
!              idof2: nominal DOF of the SPC/MPC
!
                 if(jdof1.eq.0) then
                    idof1=jdof2
                    idof2=(node1-1)*3+k
                 else
                    idof1=jdof1
                    idof2=(node2-1)*3+m
                 endif
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof2,nmpc,id)
                    if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
!
!                    regular DOF / MPC
!
                       id=ilmpc(id)
                       ist=ipompc(id)
                       index=nodempc(3,ist)
		       if(index.eq.0) cycle
                       do
                          idof2=nactdof(nodempc(2,index),
     &                                  nodempc(1,index))
                          if(idof2.ne.0) then
                             call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                          endif
                          index=nodempc(3,index)
                          if(index.eq.0) exit
                       enddo
                       cycle
                    endif
                 endif
!
              else
                 idof1=(node1-1)*3+k
                 idof2=(node2-1)*3+m
                 mpc1=0
                 mpc2=0
                 if(nmpc.gt.0) then
                    call nident(ikmpc,idof1,nmpc,id1)
                    if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
                    call nident(ikmpc,idof2,nmpc,id2)
                    if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
                 endif
                 if((mpc1.eq.1).and.(mpc2.eq.1)) then
                    id1=ilmpc(id1)
                    id2=ilmpc(id2)
                    if(id1.eq.id2) then
!
!                    MPC id1 / MPC id1
!
                       ist=ipompc(id1)
                       index1=nodempc(3,ist)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          index2=index1
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                         call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                     else
!
!                    MPC id1 / MPC id2
!
                       ist1=ipompc(id1)
                       index1=nodempc(3,ist1)
		       if(index1.eq.0) cycle
                       do
                          idof1=nactdof(nodempc(2,index1),
     &                                nodempc(1,index1))
                          ist2=ipompc(id2)
                          index2=nodempc(3,ist2)
                          if(index2.eq.0) then
                             index1=nodempc(3,index1)
                             if(index1.eq.0) then
                                exit
                             else
                                cycle
                             endif
                          endif
                          do
                             idof2=nactdof(nodempc(2,index2),
     &                                   nodempc(1,index2))
                             if((idof1.ne.0).and.(idof2.ne.0))
     &                          call insert(ipointer,mast1,irow,
     &                                   idof1,idof2,ifree,nzs_)
                             index2=nodempc(3,index2)
                             if(index2.eq.0) exit
                          enddo
                          index1=nodempc(3,index1)
                          if(index1.eq.0) exit
                       enddo
                    endif
                 endif
              endif
            enddo
          enddo
        enddo
!
!       ordering the nonzero nodes in the SUPERdiagonal columns
!       mast1 contains the row numbers column per column,
!       irow the column numbers
!
        do i=1,neq
          itot=0
          if(ipointer(i).eq.0) then
            write(*,*) 'error in mastruct: zero column'
            stop
          endif
          istart=ipointer(i)
          do
            itot=itot+1
            ikcol(itot)=mast1(istart)
            istart=irow(istart)
            if(istart.eq.0) exit
          enddo
          call isortii(ikcol,icol,itot,kflag)
          istart=ipointer(i)
          do j=1,itot-1
            mast1(istart)=ikcol(j)
            istartold=istart
            istart=irow(istart)
            irow(istartold)=i
          enddo
          mast1(istart)=ikcol(itot)
          irow(istart)=i
        enddo
!
!       defining icol and jq
!
        nsky=0
	nsky_exp=0
        do i=2,neq
          nsky_inc=i-mast1(ipointer(i))
          if(2147483647-nsky.lt.nsky_inc) then
	     nsky_exp=nsky_exp+1
	     nsky_inc=nsky_inc-2147483647
	  endif
	  nsky=nsky+nsky_inc
        enddo
!
        write(*,*) 'number of equations'
        write(*,*) neq
        write(*,*) 'number of nonzero matrix elements'
        write(*,*) ifree
        write(*,*) 'total length of the skyline'
        write(*,*) nsky_exp,'*2147483647+',nsky
        write(*,*) 'percentage of nonzero skyline elements'
        write(*,*) real(ifree)/
     &     (real(nsky+neq)+nsky_exp*real(2147483647))*100.
        write(*,*)
!
!       new meaning of icol,j1,mast1,irow:
!       - irow is going to contain the row numbers of the SUBdiagonal
!         nonzero's, column per column
!       - mast1 contains the column numbers
!       - icol(i)=# SUBdiagonal nonzero's in column i
!       - jq(i)= location in field irow of the first SUBdiagonal
!         nonzero in column i
!
        nmast=ifree
!
!       switching from a SUPERdiagonal inventary to a SUBdiagonal one
!
        call isortii(mast1,irow,nmast,kflag)
!
!       filtering out the diagonal elements and generating icol and jq
!
        isubtract=0
        do i=1,neq
           icol(i)=0
        enddo
        k=0
        do i=1,nmast
           if(mast1(i).eq.irow(i)) then
              isubtract=isubtract+1
           else
              mast1(i-isubtract)=mast1(i)
              irow(i-isubtract)=irow(i)
              if(k.ne.mast1(i)) then
                 do l=k+1,mast1(i)
                    jq(l)=i-isubtract
                 enddo
                 k=mast1(i)
              endif
              icol(k)=icol(k)+1
           endif
        enddo
        nmast=nmast-isubtract
        do l=k+1,neq+1
           jq(l)=nmast+1
        enddo
!
        do i=1,neq
           if(jq(i+1)-jq(i).gt.0) then
              call isortii(irow(jq(i)),mast1(jq(i)),jq(i+1)-jq(i),
     &          kflag)

           endif
        enddo
!
        nzs=jq(neq)-1
!
      endif
!
      return
      end

      */

invivomastruct_mastruct_705

Real time: 4.270208e+06 -- Predicted time: 4.277338e+06

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	4.277338e+06	4.278852e+06	0.0

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivoSPOOLES_Utilities_src_Utilities_newsort_IVqsortUp_975

Real time: 1.825840e+06 -- Predicted time: 1.007780e+07

Callcount: 12588

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
351	### 1	3448.46031746032	9.540000e+02	2.520000e+02	73.5
33	### 2	9311.34482758621	7.290000e+02	1.160000e+02	84.0

                *DECK ISORT
      SUBROUTINE ISORTII (IX, IY, N, KFLAG)
C***BEGIN PROLOGUE  ISORT
C***PURPOSE  Sort an array and optionally make the same interchanges in
C            an auxiliary array.  The array may be sorted in increasing
C            or decreasing order.  A slightly modified QUICKSORT
C            algorithm is used.
C***LIBRARY   SLATEC
C***CATEGORY  N6A2A
C***TYPE      INTEGER (SSORT-S, DSORT-D, ISORT-I)
C***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
C***AUTHOR  Jones, R. E., (SNLA)
C           Kahaner, D. K., (NBS)
C           Wisniewski, J. A., (SNLA)
C***DESCRIPTION
C
C   ISORT sorts array IX and optionally makes the same interchanges in
C   array IY.  The array IX may be sorted in increasing order or
C   decreasing order.  A slightly modified quicksort algorithm is used.
C
C   Description of Parameters
C      IX - integer array of values to be sorted
C      IY - integer array to be (optionally) carried along
C      N  - number of values in integer array IX to be sorted
C      KFLAG - control parameter
C            =  2  means sort IX in increasing order and carry IY along.
C            =  1  means sort IX in increasing order (ignoring IY)
C            = -1  means sort IX in decreasing order (ignoring IY)
C            = -2  means sort IX in decreasing order and carry IY along.
C
C***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
C                 for sorting with minimal storage, Communications of
C                 the ACM, 12, 3 (1969), pp. 185-187.
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   761118  DATE WRITTEN
C   810801  Modified by David K. Kahaner.
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891009  Removed unreferenced statement labels.  (WRB)
C   891009  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   901012  Declared all variables; changed X,Y to IX,IY. (M. McClain)
C   920501  Reformatted the REFERENCES section.  (DWL, WRB)
C   920519  Clarified error messages.  (DWL)
C   920801  Declarations section rebuilt and code restructured to use
C           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
C***END PROLOGUE  ISORT
C     .. Scalar Arguments ..
      INTEGER KFLAG, N
C     .. Array Arguments ..
      INTEGER IX(*), IY(*)
C     .. Local Scalars ..
      REAL R
      INTEGER I, IJ, J, K, KK, L, M, NN, T, TT, TTY, TY
C     .. Local Arrays ..
      INTEGER IL(21), IU(21)
C     .. External Subroutines ..
!      EXTERNAL XERMSG
C     .. Intrinsic Functions ..
      INTRINSIC ABS, INT
C***FIRST EXECUTABLE STATEMENT  ISORT
      NN = N
      IF (NN .LT. 1) THEN
!         CALL XERMSG ('SLATEC', 'ISORT',
!     +      'The number of values to be sorted is not positive.', 1, 1)
         RETURN
      ENDIF
C
      KK = ABS(KFLAG)
      IF (KK.NE.1 .AND. KK.NE.2) THEN
!         CALL XERMSG ('SLATEC', 'ISORT',
!     +      'The sort control parameter, K, is not 2, 1, -1, or -2.', 2,
!     +      1)
         RETURN
      ENDIF
C
C     Alter array IX to get decreasing order if needed
C
      IF (KFLAGS .LE. -1) THEN
         DO 10 I=1,NN
            IX(I) = -IX(I)
   10    CONTINUE
      ENDIF
C
      IF (KK .EQ. 2) GO TO 100
C
C     Sort IX only
C
      M = 1
      I = 1
      J = NN
      R = 0.375E0
C
   20 IF (I .EQ. J) GO TO 60
      IF (R .LE. 0.5898437E0) THEN
         R = R+3.90625E-2
      ELSE
         R = R-0.21875E0
      ENDIF
C
   30 K = I
C
C     Select a central element of the array and save it in location T
C
      IJ = I + INT((J-I)*R)
      T = IX(IJ)
C
C     If first element of array is greater than T, interchange with T
C
      IF (IX(I) .GT. T) THEN
         IX(IJ) = IX(I)
         IX(I) = T
         T = IX(IJ)
      ENDIF
      L = J
C
C     If last element of array is less than than T, interchange with T
C
      IF (IX(J) .LT. T) THEN
         IX(IJ) = IX(J)
         IX(J) = T
         T = IX(IJ)
C
C        If first element of array is greater than T, interchange with T
C
         IF (IX(I) .GT. T) THEN
            IX(IJ) = IX(I)
            IX(I) = T
            T = IX(IJ)
         ENDIF
      ENDIF
C
C     Find an element in the second half of the array which is smaller
C     than T
C
   40 L = L-1
      IF (IX(L) .GT. T) GO TO 40
C
C     Find an element in the first half of the array which is greater
C     than T
C
   50 K = K+1
      IF (IX(K) .LT. T) GO TO 50
C
C     Interchange these elements
C
      IF (K .LE. L) THEN
         TT = IX(L)
         IX(L) = IX(K)
         IX(K) = TT
         GO TO 40
      ENDIF
C
C     Save upper and lower subscripts of the array yet to be sorted
C
      IF (L-I .GT. J-K) THEN
         IL(M) = I
         IU(M) = L
         I = K
         M = M+1
      ELSE
         IL(M) = K
         IU(M) = J
         J = L
         M = M+1
      ENDIF
      GO TO 70
C
C     Begin again on another portion of the unsorted array
C
   60 M = M-1
      IF (M .EQ. 0) GO TO 190
      I = IL(M)
      J = IU(M)
C
   70 IF (J-I .GE. 1) GO TO 30
      IF (I .EQ. 1) GO TO 20
      I = I-1
C
   80 I = I+1
      IF (I .EQ. J) GO TO 60
      T = IX(I+1)
      IF (IX(I) .LE. T) GO TO 80
      K = I
C
   90 IX(K+1) = IX(K)
      K = K-1
      IF (T .LT. IX(K)) GO TO 90
      IX(K+1) = T
      GO TO 80
C
C     Sort IX and carry IY along
C
  100 M = 1
      I = 1
      J = NN
      R = 0.375E0
C
  110 IF (I .EQ. J) GO TO 150
      IF (R .LE. 0.5898437E0) THEN
         R = R+3.90625E-2
      ELSE
         R = R-0.21875E0
      ENDIF
C
  120 K = I
C
C     Select a central element of the array and save it in location T
C
      IJ = I + INT((J-I)*R)
      T = IX(IJ)
      TY = IY(IJ)
C
C     If first element of array is greater than T, interchange with T
C
      IF (IX(I) .GT. T) THEN
         IX(IJ) = IX(I)
         IX(I) = T
         T = IX(IJ)
         IY(IJ) = IY(I)
         IY(I) = TY
         TY = IY(IJ)
      ENDIF
      L = J
C
C     If last element of array is less than T, interchange with T
C
      IF (IX(J) .LT. T) THEN
         IX(IJ) = IX(J)
         IX(J) = T
         T = IX(IJ)
         IY(IJ) = IY(J)
         IY(J) = TY
         TY = IY(IJ)
C
C        If first element of array is greater than T, interchange with T
C
         IF (IX(I) .GT. T) THEN
            IX(IJ) = IX(I)
            IX(I) = T
            T = IX(IJ)
            IY(IJ) = IY(I)
            IY(I) = TY
            TY = IY(IJ)
         ENDIF
      ENDIF
C
C     Find an element in the second half of the array which is smaller
C     than T
C
  130 L = L-1
      IF (IX(L) .GT. T) GO TO 130
C
C     Find an element in the first half of the array which is greater
C     than T
C
  140 K = K+1
      IF (IX(K) .LT. T) GO TO 140
C
C     Interchange these elements
C
      IF (K .LE. L) THEN
         TT = IX(L)
         IX(L) = IX(K)
         IX(K) = TT
         TTY = IY(L)
         IY(L) = IY(K)
         IY(K) = TTY
         GO TO 130
      ENDIF
C
C     Save upper and lower subscripts of the array yet to be sorted
C
      IF (L-I .GT. J-K) THEN
         IL(M) = I
         IU(M) = L
         I = K
         M = M+1
      ELSE
         IL(M) = K
         IU(M) = J
         J = L
         M = M+1
      ENDIF
      GO TO 160
C
C     Begin again on another portion of the unsorted array
C
  150 M = M-1
      IF (M .EQ. 0) GO TO 190
      I = IL(M)
      J = IU(M)
C
  160 IF (J-I .GE. 1) GO TO 120
      IF (I .EQ. 1) GO TO 110
      I = I-1
C
  170 I = I+1
      IF (I .EQ. J) GO TO 150
      T = IX(I+1)
      TY = IY(I+1)
      IF (IX(I) .LE. T) GO TO 170
      K = I
C
  180 IX(K+1) = IX(K)
      IY(K+1) = IY(K)
      K = K-1
      IF (T .LT. IX(K)) GO TO 180
      IX(K+1) = T
      IY(K+1) = TY
      GO TO 170
C
C     Clean up
C
  190 IF (KFLAG .LE. -1) THEN
         DO 200 I=1,NN
            IX(I) = -IX(I)
  200    CONTINUE
      ENDIF
      RETURN
      END

invivoisortii_isortii__199

Real time: 1.202479e+07 -- Predicted time: 9.773259e+06

Callcount: 1441

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
474	### 1	24.5491740299654	1.421200e+04	2.082400e+04	31.7
16	### 2	215.440677966102	8.370000e+02	4.720000e+02	43.6
160	### 3	149.790322580645	2.488000e+03	2.480000e+03	0.3
175	### 4	208.017130620985	3.358000e+03	3.736000e+03	10.1
70	### 5	273.447891805887	4.447000e+03	5.028000e+03	11.5
774	### 6	186.255140186916	5.339000e+03	8.560000e+03	37.6
968	### 7	305.896466212027	2.820000e+03	6.452000e+03	56.2
285	### 8	138.059778597786	8.688000e+03	1.084000e+04	19.8
722	### 9	1	3.900312e+06	4.022132e+06	3.0

                !
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nelemprint,
     &  neprint,stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,noelplab,nodeflab,eei,een,iperturb,f,fn,
     &  nactdof,iout,qa,noprint,nodeprint,vold,b,nodeboun,ndirboun,
     &  xboun,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,vmax,neq,
     &  veold,accold,bet,gam,dtime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstiff,xstate,npmat_,epn,matname,mint_,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,sti,
     &  xstaten,eeiini,enerini,cocon,ncocon)
!
!     calculates and prints the displacements, temperatures and forces 
!     at the nodes and the stress and  strain  at the reduced integration 
!     points and at the nodes
!
!     iout=-1: v is assumed to be known and is used to
!              calculate strains, stresses..., no result output
!     iout=0: v is calculated from the system solution
!             and strains, stresses.. are calculated, no result output
!     iout=1:  v is calculated from the system solution and strains,
!              stresses.. are calculated, requested results output
!     iout=2: v is assumed to be known and is used to 
!             calculate strains, stresses..., requested results output
!
      implicit none
!
      logical calcul_fn,calcul_f,calcul_cauchy,calcul_qa,cauchy
!
      character*4 noelplab(*),nodeflab(*)
      character*8 lakon(*)
      character*20 amat,matname(*),labmpc(*)
!
      integer kon(*),konl(20),inum(*),nelemprint(*),
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(*),ielorien(*),
     &  ntmat_,ipkon(*),nactdof(0:3,*),nodeprint(*),nodeboun(*),
     &  ndirboun(*),ipompc(*),nodempc(3,*),nelem,ikboun(*),ilboun(*),
     &  ncocon(2,*)
!
      integer nk,ne,neprint,mattyp,ithermal,iprestr,i,j,k,m1,m2,jj,
     &  i1,m3,m4,kk,iener,indexe,nope,norien,iperturb,iout,noprint,
     &  nal,icmd,ihyper,nboun,nmpc,nmethod,ist,ndir,node,index,
     &  neq,kode,imat,mint3d,mint_,nopev,nfield,ndim,iorien,ielas,
     &  istiff,icmdl,ncmat_,nstate_,incrementalmpc,lb,l,jmin,jmax
!
      integer nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_
!
      real*8 co(3,*),v(0:3,*),shp(4,20),stiini(6,mint_,*),
     &  stx(6,mint_,*),stn(6,*),xl(3,20),vl(0:3,20),stre(9),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),
     &  alcon(0:6,ntmat_,*),vini(0:3,*),
     &  alzero(*),orab(7,*),elas(21),alph(6),rho,f(*),fn(0:3,*),
     &  skl(3,3),anisox(3,3,3,3),beta(6),q(0:3,20),vkl(0:3,3),
     &  t0(*),t1(*),prestr(6,*),eei(6,mint_,*),een(6,*),ckl(3,3),
     &  vold(0:3,*),b(*),xboun(*),coefmpc(*),vmax,eloc(9),veold(0:3,*),
     &  accold(0:3,*),elconloc(21),eth(6),xkl(3,3),voldl(0:3,20),epn(*),
     &  xikl(3,3),ener(mint_,*),enern(*),sti(6,mint_,*),
     &  eeiini(6,mint_,*),enerini(mint_,*),cocon(0:6,ntmat_,*)
!
      real*8 e,un,al,um,am1,am2,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
     &  xsj,qa,vj,t0l,t1l,bet,gam,dtime,forcempc,scal1,scal2,bnac,
     &  fixed_disp,weight,pgauss(3),vij,xkapp(6)
!
      real*8 plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(21,mint_,*),xstate(nstate_,mint_,*),plconloc(82),
     &  vokl(3,3),xstateini(nstate_,mint_,*),vikl(3,3),
     &  xstaten(nstate_,*)
!
      real*8 gauss2d1(2,1),gauss2d2(2,4),gauss2d3(2,9),gauss2d4(2,1),
     &  gauss2d5(2,3),gauss3d1(3,1),gauss3d2(3,8),gauss3d3(3,27),
     &  gauss3d4(3,1),gauss3d5(3,4),gauss3d6(3,15),gauss3d7(3,2),
     &  gauss3d8(3,9),gauss3d9(3,18),weight2d1(1),weight2d2(4),
     &  weight2d3(9),weight2d4(1),weight2d5(3),weight3d1(1),
     &  weight3d2(8),weight3d3(27),weight3d4(1),weight3d5(4),
     &  weight3d6(15),weight3d7(2),weight3d8(9),weight3d9(18)
!
      include "gauss.f"
!
      if(ithermal.le.1) then
         jmin=1
         jmax=3
      elseif(ithermal.eq.2) then
         jmin=0
         jmax=0
      else
         jmin=0
         jmax=3
      endif
!
      if((iout.ne.2).and.(iout.ne.-1)) then
!
         if((nmethod.ne.4).or.(iperturb.le.1)) then
!
!     extracting the displacement information from the solution
!
            do i=1,nk
               do j=jmin,jmax
                  if(nactdof(j,i).ne.0) then
                     v(j,i)=b(nactdof(j,i))
                  else
                     v(j,i)=0.d0
                  endif
               enddo
            enddo
!
!     for static perturbation steps v represents the incremental
!     displacements. For the total displacement vold must be added.
!
            if((iperturb.ne.0).and.(nmethod.eq.1)) then
               vmax=0.d0
               do i=1,nk
                  do j=jmin,jmax
                     vmax=max(vmax,dabs(v(j,i)))
c                     write(*,*) i,j,vold(j,i),v(j,i),vmax
                     v(j,i)=v(j,i)+vold(j,i)
                  enddo
               enddo
            endif
!
         else
!
!           direct integration dynamic step
!           b contains the acceleration increment
!
            if((ithermal.le.1).or.(ithermal.eq.3)) then
               scal1=bet*dtime*dtime
               scal2=gam*dtime
               vmax=0.d0
               do i=1,nk
                  do j=1,3
                     if(nactdof(j,i).ne.0) then
                        bnac=b(nactdof(j,i))
                        v(j,i)=vold(j,i)+scal1*bnac
                        vmax=max(vmax,dabs(scal1*bnac))
                        veold(j,i)=veold(j,i)+scal2*bnac
                        accold(j,i)=accold(j,i)+bnac
                     endif
                  enddo
               enddo
            endif
            if(ithermal.ge.2) then
               do i=1,nk
                  if(nactdof(0,i).ne.0) then
                     bnac=b(nactdof(0,i))
                     v(0,i)=vold(0,i)+bnac
                     vmax=max(vmax,dabs(bnac))
                     veold(0,i)=veold(0,i)+bnac/dtime
                  endif
               enddo
            endif
         endif
!
      endif
!
!        initialization
!
      calcul_fn=.false.
      calcul_f=.false.
      calcul_qa=.false.
      calcul_cauchy=.false.
!

!     determining which quantities have to be calculated
!
      if((iperturb.ge.2).or.((iperturb.eq.0).and.(iout.lt.0))) then
         if(iout.lt.1) then
            calcul_fn=.true.
            calcul_f=.true.
            calcul_qa=.true.
         else    
            calcul_cauchy=.true.
         endif
      endif
      if(((noelplab(5).eq.'RF').or.(nodeflab(5).eq.'RF'))
     &     .and.(iout.gt.0)) calcul_fn=.true.
!
!     initializing fn
!
      if(calcul_fn) then
         do i=1,nk
            do j=1,3
               fn(j,i)=0.d0
            enddo
         enddo
      endif
!
!     initializing f
!
      if(calcul_f) then
         do i=1,neq
            f(i)=0.d0
         enddo
      endif
!
!     inserting the boundary conditions
!
      do i=1,nboun
         if(ndirboun(i).gt.3) cycle
         fixed_disp=xboun(i)
         if((nmethod.eq.4).and.(iperturb.gt.1)) then
            ndir=ndirboun(i)
            node=nodeboun(i)
            accold(ndir,node)=(xboun(i)-v(ndir,node))/
     &        (bet*dtime*dtime)
            veold(ndir,node)=veold(ndir,node)+
     &        gam*dtime*accold(ndir,node)
         endif
         v(ndirboun(i),nodeboun(i))=fixed_disp
c         write(*,*) 'boun ',nodeboun(i),ndirboun(i),fixed_disp
      enddo
!
!     inserting the mpc information
!     the parameter incrementalmpc indicates whether the
!     incremental displacements enter the mpc or the total 
!     displacements (incrementalmpc=0)
!
      do i=1,nmpc
         if((labmpc(i)(1:20).eq.'                    ').or.
     &      (labmpc(i)(1:6).eq.'CYCLIC').or.
     &      (labmpc(i)(1:5).eq.'RIGID').or.
     &      (labmpc(i)(1:9).eq.'SUBCYCLIC')) then
            incrementalmpc=0
         else
            incrementalmpc=1
         endif
         ist=ipompc(i)
         node=nodempc(1,ist)
         ndir=nodempc(2,ist)
         index=nodempc(3,ist)
         fixed_disp=0.d0
         if(index.ne.0) then
            do
               if(incrementalmpc.eq.0) then
                  fixed_disp=fixed_disp-coefmpc(index)*
     &                 v(nodempc(2,index),nodempc(1,index))
               else
c                  write(*,102) 
c     &             nodempc(1,index),nodempc(2,index),
c     &             v(nodempc(2,index),nodempc(1,index)),
c     &             vold(nodempc(2,index),nodempc(1,index))
c 102              format('results,node,dir,newvalue,oldvalue',2i5,
c     &                2(1x,e11.4))
                  fixed_disp=fixed_disp-coefmpc(index)*
     &                 (v(nodempc(2,index),nodempc(1,index))-
     &                  vold(nodempc(2,index),nodempc(1,index)))
c                  write(*,*) fixed_disp,
c     & vold(nodempc(2,index),nodempc(1,index)),
c     & v(nodempc(2,index),nodempc(1,index))
               endif
               index=nodempc(3,index)
               if(index.eq.0) exit
            enddo
         endif
         fixed_disp=fixed_disp/coefmpc(ist)
c                  write(*,*) fixed_disp
         if(incrementalmpc.eq.1) then
            fixed_disp=fixed_disp+vold(ndir,node)
c            write(*,*) 'results ',i,vold(ndir,node),fixed_disp
         endif
         if((nmethod.eq.4).and.(iperturb.gt.1)) then
            accold(ndir,node)=(fixed_disp-v(ndir,node))/
     &        (bet*dtime*dtime)
            veold(ndir,node)=veold(ndir,node)+
     &        gam*dtime*accold(ndir,node)
         endif
         v(ndir,node)=fixed_disp
c         write(*,102) node,ndir,v(ndir,node),vold(ndir,node)
      enddo
!
!     check whether there are any strain output requests
!
      if ((noelplab(7).eq.'ENER').or.(nodeflab(7).eq.'ENER')) then
         iener=1
      else
         iener=0
      endif
!
      qa=0.d0
      nal=0
!
!     calculation of the stresses in the integration points
! 
      if((ithermal.le.1).or.(ithermal.eq.3)) then
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         imat=ielmat(i)
         amat=matname(imat)
         if(norien.gt.0) then
            iorien=ielorien(i)
         else
            iorien=0
         endif
!
c         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         if(lakon(i)(4:4).eq.'2') then
            nope=20
         elseif(lakon(i)(4:4).eq.'8') then
            nope=8
         elseif(lakon(i)(4:5).eq.'10') then
            nope=10
         elseif(lakon(i)(4:4).eq.'4') then
            nope=4
         elseif(lakon(i)(4:5).eq.'15') then
            nope=15
         else
            nope=6
         endif
!
         if(lakon(i)(4:5).eq.'8R') then
            mint3d=1
         elseif((lakon(i)(4:4).eq.'8').or.
     &          (lakon(i)(4:6).eq.'20R')) then
            mint3d=8
         elseif(lakon(i)(4:4).eq.'2') then
            mint3d=27
         elseif(lakon(i)(4:5).eq.'10') then
            mint3d=4
         elseif(lakon(i)(4:4).eq.'4') then
            mint3d=1
         elseif(lakon(i)(4:5).eq.'15') then
            mint3d=9
         else
            mint3d=2
         endif
!
         do j=1,nope
            konl(j)=kon(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl(j))
               vl(k,j)=v(k,konl(j))
               voldl(k,j)=vold(k,konl(j))
            enddo
         enddo
!
!        check for hyperelastic material
!
         if(nelcon(1,imat).lt.0) then
            ihyper=1
         else
            ihyper=0
         endif
!
!        q contains the nodal forces per element; initialisation of q
!
         if((iperturb.ge.2).or.((iperturb.eq.0).and.(iout.lt.1))) then
            do m1=1,nope
               do m2=1,3
                  q(m2,m1)=fn(m2,konl(m1))
               enddo
            enddo
         endif
!
         do jj=1,mint3d
            if(lakon(i)(4:5).eq.'8R') then
               xi=gauss3d1(1,jj)
               et=gauss3d1(2,jj)
               ze=gauss3d1(3,jj)
               weight=weight3d1(jj)
            elseif((lakon(i)(4:4).eq.'8').or.
     &             (lakon(i)(4:6).eq.'20R'))
     &        then
               xi=gauss3d2(1,jj)
c               if(nope.eq.20) xi=xi+1.d0
               et=gauss3d2(2,jj)
               ze=gauss3d2(3,jj)
               weight=weight3d2(jj)
            elseif(lakon(i)(4:4).eq.'2') then
c               xi=gauss3d3(1,jj)+1.d0
               xi=gauss3d3(1,jj)
               et=gauss3d3(2,jj)
               ze=gauss3d3(3,jj)
               weight=weight3d3(jj)
            elseif(lakon(i)(4:5).eq.'10') then
               xi=gauss3d5(1,jj)
               et=gauss3d5(2,jj)
               ze=gauss3d5(3,jj)
               weight=weight3d5(jj)
            elseif(lakon(i)(4:4).eq.'4') then
               xi=gauss3d4(1,jj)
               et=gauss3d4(2,jj)
               ze=gauss3d4(3,jj)
               weight=weight3d4(jj)
            elseif(lakon(i)(4:5).eq.'15') then
               xi=gauss3d8(1,jj)
               et=gauss3d8(2,jj)
               ze=gauss3d8(3,jj)
               weight=weight3d8(jj)
            else
               xi=gauss3d7(1,jj)
               et=gauss3d7(2,jj)
               ze=gauss3d7(3,jj)
               weight=weight3d7(jj)
            endif
!
            if(nope.eq.20) then
               call shape20h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp)
            else
               call shape6w(xi,et,ze,xl,xsj,shp)
            endif
!
!                 vkl(m2,m3) contains the derivative of the m2-
!                 component of the displacement with respect to
!                 direction m3
!
            do m2=1,3
               do m3=1,3
                  vkl(m2,m3)=0.d0
               enddo
            enddo
!
            do m1=1,nope
               do m2=1,3
                  do m3=1,3
                     vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
                  enddo
               enddo
            enddo
!
            kode=nelcon(1,imat)
!
!                 calculating the strain
!
!                 attention! exy,exz and eyz are engineering strains!
!
            exx=vkl(1,1)
            eyy=vkl(2,2)
            ezz=vkl(3,3)
            exy=vkl(1,2)+vkl(2,1)
            exz=vkl(1,3)+vkl(3,1)
            eyz=vkl(2,3)+vkl(3,2)
!
!                 Lagrangian/Eulerian strain tensor for geometric
!                 nonlinear calculations
!
            if(iperturb.ge.2) then
!
!                    Lagrangian strain
!
               exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
               eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
               ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
               exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
     &              vkl(3,1)*vkl(3,2)
               exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
     &              vkl(3,1)*vkl(3,3)
               eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
     &              vkl(3,2)*vkl(3,3)
!
            endif
!
!                 storing the local strains
!
            if((iperturb.ne.0).or.(iout.ge.0)) then
               eloc(1)=exx
               eloc(2)=eyy
               eloc(3)=ezz
               eloc(4)=exy/2.d0
               eloc(5)=exz/2.d0
               eloc(6)=eyz/2.d0
            else
!
!                    linear iteration within a nonlinear increment:
!                    tangent matrix is calculated at strain at
!                    the end of the previous increment
!
               do m1=1,6
                  eloc(m1)=eei(m1,jj,i)
               enddo
            endif
!
!                 calculating the deformation gradient (needed to
!                 convert the element stiffness matrix from spatial
!                 coordinates to material coordinates
!                 deformation plasticity)
!
            if((kode.eq.-50).or.(kode.le.-100)) then
!
!                    calculating the deformation gradient
!
               xkl(1,1)=vkl(1,1)+1
               xkl(2,2)=vkl(2,2)+1.
               xkl(3,3)=vkl(3,3)+1.
               xkl(1,2)=vkl(1,2)
               xkl(1,3)=vkl(1,3)
               xkl(2,3)=vkl(2,3)
               xkl(2,1)=vkl(2,1)
               xkl(3,1)=vkl(3,1)
               xkl(3,2)=vkl(3,2)
!
!                    calculating the Jacobian
!
               vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
     &              -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
     &              +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
!
!              inversion of the deformation gradient (only for
!              deformation plasticity)
!
               if(kode.eq.-50) then
!
                  ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
                  ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
                  ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
                  ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
                  ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
                  ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
                  ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
                  ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
                  ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
!
!                 converting the Lagrangian strain into Eulerian
!                 strain (only for deformation plasticity)
!
                  cauchy=.false.
                  call str2mat(eloc,ckl,vj,cauchy)
               endif
!                        
            endif
!
!                 calculating fields for incremental plasticity
!
            if(kode.le.-100) then
!
!              calculating the deformation gradient at the
!              start of the increment
!
!              calculating the displacement gradient at the
!              start of the increment
!
               do m2=1,3
                  do m3=1,3
                     vikl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vikl(m2,m3)=vikl(m2,m3)
     &                       +shp(m3,m1)*vini(m2,konl(m1))
                     enddo
                  enddo
               enddo
!
!              calculating the deformation gradient of the old
!              fields
!
               xikl(1,1)=vikl(1,1)+1
               xikl(2,2)=vikl(2,2)+1.
               xikl(3,3)=vikl(3,3)+1.
               xikl(1,2)=vikl(1,2)
               xikl(1,3)=vikl(1,3)
               xikl(2,3)=vikl(2,3)
               xikl(2,1)=vikl(2,1)
               xikl(3,1)=vikl(3,1)
               xikl(3,2)=vikl(3,2)
!
!              calculating the Jacobian
!
               vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,2))
     &              -xikl(1,2)*(xikl(2,1)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,1))
     &              +xikl(1,3)*(xikl(2,1)*xikl(3,2)
     &              -xikl(2,2)*xikl(3,1))
!
!              stresses at the start of the increment
!
               do m1=1,6
                  stre(m1)=stiini(m1,jj,i)
               enddo
!
            endif
!
!                 prestress values
!
            if(iprestr.eq.0) then
               do kk=1,6
                  beta(kk)=0.d0
               enddo
            else
               do kk=1,6
                  beta(kk)=-prestr(kk,i)
               enddo
            endif
!
            if(ithermal.ge.1) then
!
!              calculating the temperature difference in
!              the integration point
!
               t0l=0.d0
               t1l=0.d0
               if(ithermal.eq.1) then
                  if(lakon(i)(4:5).eq.'8 ') then
                     do i1=1,nope
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+t1(konl(i1))/8.d0
                     enddo
                  elseif(lakon(i)(4:6).eq.'20 ') then
                     call lintemp(t0,t1,konl,nope,jj,t0l,t1l)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*t1(konl(i1))
                     enddo
                  endif
               elseif(ithermal.ge.2) then
                  if(lakon(i)(4:5).eq.'8 ') then
                     do i1=1,nope
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+vold(0,konl(i1))/8.d0
                     enddo
                  elseif(lakon(i)(4:6).eq.'20 ') then
                     call lintemp_th(t0,vold,konl,nope,jj,t0l,t1l)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*vold(0,konl(i1))
                     enddo
                  endif
               endif
               tt=t1l-t0l
            endif
!
!                 calculating the coordinates of the integration point
!                 for material orientation purposes (for cylindrical
!                 coordinate systems)
!
            if(iorien.gt.0) then
               do j=1,3
                  pgauss(j)=0.d0
                  do i1=1,nope
                     pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
                  enddo
               enddo
            endif
!
!                 material data; for linear elastic materials
!                 this includes the calculation of the stiffness
!                 matrix
!
            istiff=0
!
            call materialdata(elcon,nelcon,rhcon,nrhcon,alcon,
     &           nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
     &           elas,alph,rho,i,ithermal,alzero,mattyp,
     &           t0l,t1l,ihyper,istiff,
     &           elconloc,eth,kode,plicon,nplicon,
     &           plkcon,nplkcon,npmat_,
     &           plconloc,mint_,dtime,i,jj,
     &           xstiff,ncmat_)
!
            if((ihyper.eq.0).and.
     &           ((iperturb.ne.0).or.(iout.ge.0))) then
               if(mattyp.eq.1) then
                  e=elas(1)
                  un=elas(2)
                  al=un*e/(1.d0+un)/(1.d0-2.d0*un)
                  um=e/2.d0/(1.d0+un)
                  am1=al+2.d0*um
                  am2=2.d0*um
               elseif(mattyp.eq.3) then
                  call anisotropic(elas,anisox)
               endif
!
               icmdl=3
               call linel(ithermal,mattyp,beta,al,um,am1,alph,tt,
     &              elas,icmdl,exx,eyy,ezz,exy,exz,eyz,stre,
     &              anisox)
!
            elseif(ihyper.ne.0) then
               mattyp=3
               call hyperel(elconloc,elas,eloc,kode,eth,ithermal,
     &              icmd,beta,stre,xkl,ckl,vj,xikl,vij,
     &              plconloc,xstate,xstateini,ielas,
     &              amat,t1l,dtime,i,jj,nstate_,mint_,
     &              iorien,pgauss,orab)
               do m1=1,21
                  xstiff(m1,jj,i)=elas(m1)
               enddo
            endif
!
            if((iperturb.eq.0).and.(iout.lt.0)) then
!
!                    if the forced displacements were changed at
!                    the start of a nonlinear step, the nodal
!                    forces due do this displacements are 
!                    calculated in a purely linear way, and
!                    the first iteration is purely linear in order
!                    to allow the displacements to redistribute
!                    in a quasi-static way (only applies to
!                    quasi-static analyses (*STATIC))
!
               do m2=1,3
                  do m3=1,3
                     vokl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vokl(m2,m3)=vokl(m2,m3)+
     &                       shp(m3,m1)*voldl(m2,m1)
                     enddo
                  enddo
               enddo
!
               eloc(1)=exx-vokl(1,1)
               eloc(2)=eyy-vokl(2,2)
               eloc(3)=ezz-vokl(3,3)
               eloc(4)=exy-(vokl(1,2)+vokl(2,1))
               eloc(5)=exz-(vokl(1,3)+vokl(3,1))
               eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
               if(mattyp.eq.1) then
                  e=elas(1)
                  un=elas(2)
                  um=e/(1.d0+un)
                  al=un*um/(1.d0-2.d0*un)
                  um=um/2.d0
                  am1=al*(eloc(1)+eloc(2)+eloc(3))
                  stre(1)=am1+2.d0*um*eloc(1)
                  stre(2)=am1+2.d0*um*eloc(2)
                  stre(3)=am1+2.d0*um*eloc(3)
                  stre(4)=um*eloc(4)
                  stre(5)=um*eloc(5)
                  stre(6)=um*eloc(6)
               elseif(mattyp.eq.2) then
                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
     &                 +eloc(3)*elas(4)
                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
     &                 +eloc(3)*elas(5)
                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
     &                 +eloc(3)*elas(6)
                  stre(4)=eloc(4)*elas(7)
                  stre(5)=eloc(5)*elas(8)
                  stre(6)=eloc(6)*elas(9)
               elseif(mattyp.eq.3) then
                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+
     &                 eloc(3)*elas(4)+eloc(4)*elas(7)+
     &                 eloc(5)*elas(11)+eloc(6)*elas(16)
                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+
     &                 eloc(3)*elas(5)+eloc(4)*elas(8)+
     &                 eloc(5)*elas(12)+eloc(6)*elas(17)
                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+
     &                 eloc(3)*elas(6)+eloc(4)*elas(9)+
     &                 eloc(5)*elas(13)+eloc(6)*elas(18)
                  stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+
     &                 eloc(3)*elas(9)+eloc(4)*elas(10)+
     &                 eloc(5)*elas(14)+eloc(6)*elas(19)
                  stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+
     &                 eloc(3)*elas(13)+eloc(4)*elas(14)+
     &                 eloc(5)*elas(15)+eloc(6)*elas(20)
                  stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+
     &                 eloc(3)*elas(18)+eloc(4)*elas(19)+
     &                 eloc(5)*elas(20)+eloc(6)*elas(21)
               endif
            endif
! 
!           updating the energy
!
            if(iener.eq.1) then
               ener(jj,i)=enerini(jj,i)+
     &           ((eloc(1)-eeiini(1,jj,i))*(stre(1)+stiini(1,jj,i))+
     &            (eloc(2)-eeiini(2,jj,i))*(stre(2)+stiini(2,jj,i))+
     &            (eloc(3)-eeiini(3,jj,i))*(stre(3)+stiini(3,jj,i)))
     &            /2.d0+
     &            (eloc(4)-eeiini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
     &            (eloc(5)-eeiini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
     &            (eloc(6)-eeiini(6,jj,i))*(stre(6)+stiini(6,jj,i))
c               write(*,*) 'energy',jj,i,ener(jj,i)
            endif
!
            eei(1,jj,i)=eloc(1)
            eei(2,jj,i)=eloc(2)
            eei(3,jj,i)=eloc(3)
            eei(4,jj,i)=eloc(4)
            eei(5,jj,i)=eloc(5)
            eei(6,jj,i)=eloc(6)
!
            skl(1,1)=stre(1)
            skl(2,2)=stre(2)
            skl(3,3)=stre(3)
            skl(2,1)=stre(4)
            skl(3,1)=stre(5)
            skl(3,2)=stre(6)
!
            stx(1,jj,i)=skl(1,1)
            stx(2,jj,i)=skl(2,2)
            stx(3,jj,i)=skl(3,3)
            stx(4,jj,i)=skl(2,1)
            stx(5,jj,i)=skl(3,1)
            stx(6,jj,i)=skl(3,2)
!
            skl(1,2)=skl(2,1)
            skl(1,3)=skl(3,1)
            skl(2,3)=skl(3,2)
!
!                 calculation of the nodal forces
!
            if(calcul_fn)then
!
!                    calculating fn using skl
!
               do m1=1,nope
                  do m2=1,3
!
!                          linear elastic part
!                           
                     do m3=1,3
                        fn(m2,konl(m1))=fn(m2,konl(m1))+
     &                       xsj*skl(m2,m3)*shp(m3,m1)*weight
                     enddo
!
!                          nonlinear geometric part
!
                     if(iperturb.ge.2) then
                        do m3=1,3
                           do m4=1,3
                              fn(m2,konl(m1))=fn(m2,konl(m1))+
     &                             xsj*skl(m4,m3)*weight*
     &                             (vkl(m2,m4)*shp(m3,m1)+
     &                             vkl(m2,m3)*shp(m4,m1))/2.d0
                           enddo
                        enddo
                     endif
!
                  enddo
               enddo
            endif
!
!           calculation of the Cauchy stresses
!
            if(calcul_cauchy) then
!
!              changing the displacement gradients into
!              deformation gradients
!
               if(kode.ne.-50) then
                  xkl(1,1)=vkl(1,1)+1
                  xkl(2,2)=vkl(2,2)+1.
                  xkl(3,3)=vkl(3,3)+1.
                  xkl(1,2)=vkl(1,2)
                  xkl(1,3)=vkl(1,3)
                  xkl(2,3)=vkl(2,3)
                  xkl(2,1)=vkl(2,1)
                  xkl(3,1)=vkl(3,1)
                  xkl(3,2)=vkl(3,2)
!
                  vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
     &                 -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
     &                 +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
               endif
!
               do m1=1,3
                  do m2=1,m1
                     ckl(m1,m2)=0.d0
                     do m3=1,3
                        do m4=1,3
                           ckl(m1,m2)=ckl(m1,m2)+
     &                          skl(m3,m4)*xkl(m1,m3)*xkl(m2,m4)
                        enddo
                     enddo
                     ckl(m1,m2)=ckl(m1,m2)/vj
                  enddo
               enddo
!
               stx(1,jj,i)=ckl(1,1)
               stx(2,jj,i)=ckl(2,2)
               stx(3,jj,i)=ckl(3,3)
               stx(4,jj,i)=ckl(2,1)
               stx(5,jj,i)=ckl(3,1)
               stx(6,jj,i)=ckl(3,2)
            endif
!
         enddo
!
!        q contains the contributions to the nodal force in the nodes
!        belonging to the element at stake from other elements (elements
!        already treated). These contributions have to be
!        subtracted to get the contributions attributable to the element
!        at stake only
!
         if(calcul_qa) then
            do m1=1,nope
               do m2=1,3
                  qa=qa+dabs(fn(m2,konl(m1))-q(m2,m1))
               enddo
            enddo
            nal=nal+3*nope
         endif
      enddo
!
      endif
!
!     calculation of temperatures and thermal flux
!
      if(ithermal.ge.2) then
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         imat=ielmat(i)
         amat=matname(imat)
         if(norien.gt.0) then
            iorien=ielorien(i)
         else
            iorien=0
         endif
!
         indexe=ipkon(i)
         if(lakon(i)(4:4).eq.'2') then
            nope=20
         elseif(lakon(i)(4:4).eq.'8') then
            nope=8
         elseif(lakon(i)(4:5).eq.'10') then
            nope=10
         elseif(lakon(i)(4:4).eq.'4') then
            nope=4
         elseif(lakon(i)(4:5).eq.'15') then
            nope=15
         else
            nope=6
         endif
!
         if(lakon(i)(4:5).eq.'8R') then
            mint3d=1
         elseif((lakon(i)(4:4).eq.'8').or.
     &          (lakon(i)(4:6).eq.'20R')) then
            mint3d=8
         elseif(lakon(i)(4:4).eq.'2') then
            mint3d=27
         elseif(lakon(i)(4:5).eq.'10') then
            mint3d=4
         elseif(lakon(i)(4:4).eq.'4') then
            mint3d=1
         elseif(lakon(i)(4:5).eq.'15') then
            mint3d=9
         else
            mint3d=2
         endif
!
         do j=1,nope
            konl(j)=kon(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl(j))
            enddo
            vl(0,j)=v(0,konl(j))
            voldl(0,j)=vold(0,konl(j))
         enddo
!
!        q contains the nodal forces per element; initialisation of q
!
         if((iperturb.ge.2).or.((iperturb.eq.0).and.(iout.lt.1))) then
            do m1=1,nope
               q(0,m1)=fn(0,konl(m1))
            enddo
         endif
!
         do jj=1,mint3d
            if(lakon(i)(4:5).eq.'8R') then
               xi=gauss3d1(1,jj)
               et=gauss3d1(2,jj)
               ze=gauss3d1(3,jj)
               weight=weight3d1(jj)
            elseif((lakon(i)(4:4).eq.'8').or.
     &             (lakon(i)(4:6).eq.'20R'))
     &        then
               xi=gauss3d2(1,jj)
               et=gauss3d2(2,jj)
               ze=gauss3d2(3,jj)
               weight=weight3d2(jj)
            elseif(lakon(i)(4:4).eq.'2') then
               xi=gauss3d3(1,jj)
               et=gauss3d3(2,jj)
               ze=gauss3d3(3,jj)
               weight=weight3d3(jj)
            elseif(lakon(i)(4:5).eq.'10') then
               xi=gauss3d5(1,jj)
               et=gauss3d5(2,jj)
               ze=gauss3d5(3,jj)
               weight=weight3d5(jj)
            elseif(lakon(i)(4:4).eq.'4') then
               xi=gauss3d4(1,jj)
               et=gauss3d4(2,jj)
               ze=gauss3d4(3,jj)
               weight=weight3d4(jj)
            elseif(lakon(i)(4:5).eq.'15') then
               xi=gauss3d8(1,jj)
               et=gauss3d8(2,jj)
               ze=gauss3d8(3,jj)
               weight=weight3d8(jj)
            else
               xi=gauss3d7(1,jj)
               et=gauss3d7(2,jj)
               ze=gauss3d7(3,jj)
               weight=weight3d7(jj)
            endif
!
            if(nope.eq.20) then
               call shape20h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp)
            else
               call shape6w(xi,et,ze,xl,xsj,shp)
            endif
!
!                 vkl(m2,m3) contains the derivative of the m2-
!                 component of the displacement with respect to
!                 direction m3
!
            do m3=1,3
               vkl(0,m3)=0.d0
            enddo
!
            do m1=1,nope
               do m3=1,3
                  vkl(0,m3)=vkl(0,m3)+shp(m3,m1)*vl(0,m1)
               enddo
            enddo
!
            kode=nelcon(1,imat)
!
!           storing the temperature gradients
!
            eloc(7)=vkl(0,1)
            eloc(8)=vkl(0,2)
            eloc(9)=vkl(0,3)
!
!              calculating the temperature difference in
!              the integration point
!
            t1l=0.d0
            if(lakon(i)(4:5).eq.'8 ') then
               do i1=1,nope
                  t1l=t1l+vold(0,konl(i1))/8.d0
               enddo
            elseif(lakon(i)(4:6).eq.'20 ') then
               call lintemp_th(t0,vold,konl,nope,jj,t0l,t1l)
            else
               do i1=1,nope
                  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
               enddo
            endif
!
!           calculating the coordinates of the integration point
!           for material orientation purposes (for cylindrical
!           coordinate systems)
!
            if(iorien.gt.0) then
               do j=1,3
                  pgauss(j)=0.d0
                  do i1=1,nope
                     pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
                  enddo
               enddo
            endif
!
!                 material data; for linear elastic materials
!                 this includes the calculation of the stiffness
!                 matrix
!
            call matdata_co(cocon,ncocon,imat,iorien,pgauss,orab,ntmat_,
     &           xkapp,mattyp,t1l)
!
            if(mattyp.eq.1) then
               do m1=1,3
                  stre(6+m1)=xkapp(1)*vkl(0,m1)
               enddo
            elseif(mattyp.eq.2) then
               stre(7)=xkapp(1)*vkl(0,1)
               stre(8)=xkapp(2)*vkl(0,2)
               stre(9)=xkapp(3)*vkl(0,3)
            else
               stre(7)=xkapp(1)*vkl(0,1)+xkapp(2)*vkl(0,2)+
     &                 xkapp(4)*vkl(0,3)
               stre(8)=xkapp(2)*vkl(0,1)+xkapp(3)*vkl(0,2)+
     &                 xkapp(5)*vkl(0,3)
               stre(9)=xkapp(4)*vkl(0,1)+xkapp(5)*vkl(0,2)+
     &                 xkapp(6)*vkl(0,3)
            endif
! 
!           updating the energy
!
c            if(iener.eq.1) then
c               ener(jj,i)=enerini(jj,i)+
c     &           ((eloc(1)-eeiini(1,jj,i))*(stre(1)+stiini(1,jj,i))+
c     &            (eloc(2)-eeiini(2,jj,i))*(stre(2)+stiini(2,jj,i))+
c     &            (eloc(3)-eeiini(3,jj,i))*(stre(3)+stiini(3,jj,i)))
c     &            /2.d0+
c     &            (eloc(4)-eeiini(4,jj,i))*(stre(4)+stiini(4,jj,i))+
c     &            (eloc(5)-eeiini(5,jj,i))*(stre(5)+stiini(5,jj,i))+
c     &            (eloc(6)-eeiini(6,jj,i))*(stre(6)+stiini(6,jj,i))
c               write(*,*) 'energy',jj,i,ener(jj,i)
c            endif
!
!           dimension of stx should be 9!!!!!!
!
c            stx(7,jj,i)=skl(1,1)
c            stx(8,jj,i)=skl(2,2)
c            stx(9,jj,i)=skl(3,3)
!
!           calculation of the nodal flux
!
            if(calcul_fn)then
!
!                    calculating fn using skl
!
               do m1=1,nope
                  do m3=1,3
                     fn(0,konl(m1))=fn(0,konl(m1))+
     &                    xsj*stre(m3)*shp(m3,m1)*weight
                  enddo
               enddo
            endif
         enddo
!
!        q contains the contributions to the nodal force in the nodes
!        belonging to the element at stake from other elements (elements
!        already treated). These contributions have to be
!        subtracted to get the contributions attributable to the element
!        at stake only
!
         if(calcul_qa) then
            do m1=1,nope
               do m2=1,3
                  qa=qa+dabs(fn(m2,konl(m1))-q(m2,m1))
               enddo
            enddo
            nal=nal+3*nope
         endif
      enddo
!
      endif
!
      if(calcul_qa) qa=qa/nal
!
!     subtracting the mpc force (for each linear mpc there is one
!     force; the actual force in a node belonging to the mpc is
!     obtained by multiplying this force with the nodal coefficient.
!     The force has to be subtracted from f, since it does not
!     appear on the rhs of the equations system
!
        do i=1,nmpc
            ist=ipompc(i)
            node=nodempc(1,ist)
            ndir=nodempc(2,ist)
            forcempc=fn(ndir,node)/coefmpc(ist)
            fn(ndir,node)=0.d0
            index=nodempc(3,ist)
            if(index.eq.0) cycle
            do
               node=nodempc(1,index)
               ndir=nodempc(2,index)
               fn(ndir,node)=fn(ndir,node)-coefmpc(index)*forcempc
               index=nodempc(3,index)
               if(index.eq.0) exit
            enddo
         enddo
!
!     calculating the system force vector
!
      if(calcul_f) then
         do i=1,nk
            do j=1,3
c               write(*,500) i,j,fn(j,i)
c 500           format(i5,',',i5,',',e11.4)
               if(nactdof(j,i).ne.0) then
                  f(nactdof(j,i))=fn(j,i)
               endif
            enddo
         enddo
!
!     subtracting the mpc force (for each linear mpc there is one
!     force; the actual force in a node belonging to the mpc is
!     obtained by multiplying this force with the nodal coefficient.
!     The force has to be subtracted from f, since it does not
!     appear on the rhs of the equations system
!
c         do i=1,nmpc
c            ist=ipompc(i)
c            node=nodempc(1,ist)
c            ndir=nodempc(2,ist)
c            forcempc=fn(ndir,node)/coefmpc(ist)
c            index=nodempc(3,ist)
c            if(index.eq.0) cycle
c            do
c               node=nodempc(1,index)
c               ndir=nodempc(2,index)
c               j=nactdof(ndir,node)
c
c to be modified for nodes which do not belong to any elements:
c  thus there is no stress and no force is calculated
c  example: rigid body motion
c
c               if(j.ne.0) then
c                  f(j)=f(j)-coefmpc(index)*forcempc
c               endif
c               index=nodempc(3,index)
c               if(index.eq.0) exit
c            enddo
c         enddo
      endif
!
!     no print requests
!
      if(iout.le.0) return
!
!     nodal printing requests
!
      write(*,*) 'Printing nodal quantities'
      write(*,*)
!
      if(noprint.ne.0) then
         if(noelplab(1).eq.'U   ') then
            write(5,*)
            write(5,*) 'displacements (vx,vy,vz)'
            write(5,*)
            do i=1,noprint
               write(5,'(i5,1p,3(1x,e11.4))') nodeprint(i),
     &              (v(j,nodeprint(i)),j=1,3)
            enddo
         endif
!
         if(noelplab(2).eq.'NT  ') then
            write(5,*)
            write(5,*) 'temperature'
            write(5,*)
            if(ithermal.le.1) then
               do i=1,noprint
                  write(5,'(i5,1x,1p,e11.4)') nodeprint(i),
     &                 t1(nodeprint(i))
               enddo
            else
               do i=1,noprint
                  write(5,'(i5,1x,1p,e11.4)') nodeprint(i),
     &                 vold(0,nodeprint(i))
               enddo
            endif
         endif
         
!
         if(noelplab(5).eq.'RF  ') then
            write(5,*)
            write(5,*) 'forces (fx,fy,fz)'
            write(5,*)
            do i=1,noprint
               write(5,'(i5,1p,3(1x,e11.4))') nodeprint(i),
     &              (fn(j,nodeprint(i)),j=1,3)
            enddo
         endif
      endif
!
!     elemental printing requests
!
      write(*,*) 'Printing element quantities'
      write(*,*)
!
      if(neprint.ne.0) then
         if(noelplab(3).eq.'S   ') then
            write(5,*)
            write(5,*) 
     &         'stresses (elem, integ.pnt.,sxx,syy,szz,sxy,sxz,syz)'
            write(5,*)
            do i=1,neprint
               nelem=nelemprint(i)
               if(ipkon(nelem).lt.0) cycle
               if(lakon(nelem)(4:5).eq.'8R') then
                  mint3d=1
               elseif((lakon(nelem)(4:4).eq.'8').or.
     &                (lakon(nelem)(4:6).eq.'20R')) then
                  mint3d=8
               elseif(lakon(nelem)(4:4).eq.'2') then
                  mint3d=27
               elseif(lakon(nelem)(4:5).eq.'10') then
                  mint3d=4
               elseif(lakon(nelem)(4:4).eq.'4') then
                  mint3d=1
               elseif(lakon(nelem)(4:5).eq.'15') then
                  mint3d=9
               else
                  mint3d=2
               endif
               do j=1,mint3d
                  write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                 (stx(k,j,nelem),k=1,6)
               enddo
            enddo
         endif
!
         if(noelplab(4).eq.'E   ') then
            write(5,*)
            write(5,*) 
     &         'strains (elem, integ.pnt.,exx,eyy,ezz,exy,exz,eyz)'
            write(5,*)
            do i=1,neprint
               nelem=nelemprint(i)
               if(ipkon(nelem).lt.0) cycle
               if(lakon(nelem)(4:5).eq.'8R') then
                  mint3d=1
               elseif((lakon(nelem)(4:4).eq.'8').or.
     &                (lakon(nelem)(4:6).eq.'20R')) then
                  mint3d=8
               elseif(lakon(nelem)(4:4).eq.'2') then
                  mint3d=27
               elseif(lakon(nelem)(4:5).eq.'10') then
                  mint3d=4
               elseif(lakon(nelem)(4:4).eq.'4') then
                  mint3d=1
               elseif(lakon(nelem)(4:5).eq.'15') then
                  mint3d=9
               else
                  mint3d=2
               endif
               do j=1,mint3d
                  write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                 (eei(k,j,nelem),k=1,6)
               enddo
            enddo
         endif
!
         if(noelplab(6).eq.'PE  ') then
            write(5,*)
            write(5,*) 
     &         'equivalent plastic strain (elem, integ.pnt.,pe)'
            write(5,*)
            do i=1,neprint
               nelem=nelemprint(i)
               if(ipkon(nelem).lt.0) cycle
               if(lakon(nelem)(4:5).eq.'8R') then
                  mint3d=1
               elseif((lakon(nelem)(4:4).eq.'8').or.
     &                (lakon(nelem)(4:6).eq.'20R')) then
                  mint3d=8
               elseif(lakon(nelem)(4:4).eq.'2') then
                  mint3d=27
               elseif(lakon(nelem)(4:5).eq.'10') then
                  mint3d=4
               elseif(lakon(nelem)(4:4).eq.'4') then
                  mint3d=1
               elseif(lakon(nelem)(4:5).eq.'15') then
                  mint3d=9
               else
                  mint3d=2
               endif
               do j=1,mint3d
                  write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                 xstate(1,j,nelem)
               enddo
            enddo
         endif
!
         if(noelplab(7).eq.'ENER') then
            write(5,*)
            write(5,*) 
     &         'energy density (elem, integ.pnt.,energy)'
            write(5,*)
            do i=1,neprint
               nelem=nelemprint(i)
               if(ipkon(nelem).lt.0) cycle
               if(lakon(nelem)(4:5).eq.'8R') then
                  mint3d=1
               elseif((lakon(nelem)(4:4).eq.'8').or.
     &                (lakon(nelem)(4:6).eq.'20R')) then
                  mint3d=8
               elseif(lakon(nelem)(4:4).eq.'2') then
                  mint3d=27
               elseif(lakon(nelem)(4:5).eq.'10') then
                  mint3d=4
               elseif(lakon(nelem)(4:4).eq.'4') then
                  mint3d=1
               elseif(lakon(nelem)(4:5).eq.'15') then
                  mint3d=9
               else
                  mint3d=2
               endif
               do j=1,mint3d
                  write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                 ener(j,nelem)
               enddo
            enddo
         endif
!
         if(noelplab(8).eq.'SDV ') then
            do l=1,(nstate_+5)/6
               lb=(l-1)*6
               write(5,*)
               if(l.eq.(nstate_+5)/6) then
                  write(5,300) (lb+k,k=1,nstate_-lb)
               else
                  write(5,300) (lb+k,k=1,6)
               endif
 300           format('internal state variables (elem, integ.pnt.)',
     &             6(i2))
               write(5,*)
               do i=1,neprint
                  nelem=nelemprint(i)
                  if(ipkon(nelem).lt.0) cycle
                  if(lakon(nelem)(4:5).eq.'8R') then
                     mint3d=1
                  elseif((lakon(nelem)(4:4).eq.'8').or.
     &                    (lakon(nelem)(4:6).eq.'20R')) then
                     mint3d=8
                  elseif(lakon(nelem)(4:4).eq.'2') then
                     mint3d=27
                  elseif(lakon(nelem)(4:5).eq.'10') then
                     mint3d=4
                  elseif(lakon(nelem)(4:4).eq.'4') then
                     mint3d=1
                  elseif(lakon(nelem)(4:5).eq.'15') then
                     mint3d=9
                  else
                     mint3d=2
                  endif
                  do j=1,mint3d
                     if(l.eq.(nstate_+5)/6) then
                        write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                       (xstate(lb+k,j,nelem),k=1,nstate_-lb)
                     else
                        write(5,'(2i5,1p,6(1x,e11.4))') nelem,j,
     &                       (xstate(lb+k,j,nelem),k=1,6)
                     endif
                  enddo
               enddo
            enddo
         endif
!
      endif
!
!     determining the stresses in the nodes for output in frd format
!
      if(nodeflab(3).eq.'S   ') then
         nfield=6
         ndim=6
         call extrapolate(stx,stn,ipkon,inum,kon,lakon,nfield,nk,
     &        ne,mint_,ndim)
      endif
!
!     determining the strains in the nodes for output in frd format
!
      if(nodeflab(4).eq.'E   ') then
         nfield=6
         ndim=6
         call extrapolate(eei,een,ipkon,inum,kon,lakon,nfield,nk,
     &        ne,mint_,ndim)
      endif
!
!     determining the plastic equivalent strain in the nodes 
!     for output in frd format
!
      if(nodeflab(6).eq.'PE  ') then
         nfield=1
         ndim=nstate_
         call extrapolate(xstate,epn,ipkon,inum,kon,lakon,nfield,nk,
     &        ne,mint_,ndim)
      endif
!
!     determining the total energy in the nodes 
!     for output in frd format
!
      if(nodeflab(7).eq.'ENER') then
         nfield=1
         ndim=1
         call extrapolate(ener,enern,ipkon,inum,kon,lakon,nfield,nk,
     &        ne,mint_,ndim)
      endif
!
!     determining the internal state variables in the nodes 
!     for output in frd format
!
      if(nodeflab(8).eq.'SDV ') then
         nfield=nstate_
         ndim=nstate_
         call extrapolate(xstate,xstaten,ipkon,inum,kon,lakon,nfield,nk,
     &        ne,mint_,ndim)
      endif
!
!     if neither stresses nor strains are requested: calculate
!     inum: used in subroutine frd to determine which nodes are
!     active in the model (there is always frd output, e.g. U or
!     NT; default is U and S)
!
      if((nodeflab(3).ne.'S   ').and.(nodeflab(4).ne.'E   ').and.
     &   (nodeflab(6).ne.'PE  ').and.(nodeflab(7).ne.'ENER').and.
     &   (nodeflab(8).ne.'SDV ')) then
!
         do i=1,nk
            inum(i)=0
         enddo
!
         do i=1,ne
!
            if(ipkon(i).lt.0) cycle
            indexe=ipkon(i)
            if(lakon(i)(1:4).eq.'C3D1') then
               nopev=4
            else
               nopev=8
            endif
!
            do j=1,nopev
               inum(kon(indexe+j))=inum(kon(indexe+j))+1
            enddo
         enddo
!
!        determining the midside nodes
!
         do i=1,ne
!
            if(ipkon(i).lt.0) cycle
            indexe=ipkon(i)
            if(lakon(i)(4:4).eq.'2') then
               nope=20
               nopev=8
            elseif(lakon(i)(4:4).eq.'8') then
               nope=8
               nopev=8
            elseif(lakon(i)(4:5).eq.'10') then
               nope=10
               nopev=4
            elseif(lakon(i)(4:4).eq.'4') then
               nope=4
               nopev=4
            elseif(lakon(i)(4:5).eq.'15') then
               nope=15
               nopev=6
            else
               nope=6
               nopev=6
            endif
!
            do j=nopev,nope
               if(inum(kon(indexe+j)).eq.0) inum(kon(indexe+j))=-1
            enddo
         enddo
      endif
!
      return
      end

invivoresults_results__1301

Real time: 4.682592e+06 -- Predicted time: 1.284400e+05

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	1.284400e+05	4.781560e+06	97.3

                CODE NOT FOUND -> THIS CODELET NOT IN regions.csv?

invivoSPOOLES_Utilities_src_Utilities_sortAndCompress_IV2DVsortUpAndCompress_339

Real time: 3.298192e+06 -- Predicted time: 3.343572e+06

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	3.343572e+06	3.358316e+06	0.4

                !
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!     
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of 
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nelemprint,
     &  neprint,stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
     &  iprestr,noelplab,nodeflab,eei,een,iperturb,f,fn,
     &  nactdof,iout,qa,noprint,nodeprint,vold,b,nodeboun,ndirboun,
     &  xboun,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,vmax,neq,
     &  veold,accold,bet,gam,dtime,plicon,nplicon,plkcon,nplkcon,
     &  xstateini,xstiff,xstate,npmat_,epn,matname,mint_,ielas,icmd,
     &  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,sti,
     &  xstaten,eeiini,enerini,cocon,ncocon)
!
!     calculates and prints the displacements, temperatures and forces 
!     at the nodes and the stress and  strain  at the reduced integration 
!     points and at the nodes
!
!     iout=-1: v is assumed to be known and is used to
!              calculate strains, stresses..., no result output
!     iout=0: v is calculated from the system solution
!             and strains, stresses.. are calculated, no result output
!     iout=1:  v is calculated from the system solution and strains,
!              stresses.. are calculated, requested results output
!     iout=2: v is assumed to be known and is used to 
!             calculate strains, stresses..., requested results output
!
      implicit none
!
      logical calcul_fn,calcul_f,calcul_cauchy,calcul_qa,cauchy
!
      character*4 noelplab(*),nodeflab(*)
      character*8 lakon(*)
      character*20 amat,matname(*),labmpc(*)
!
      integer kon(*),konl(20),inum(*),nelemprint(*),
     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(*),ielorien(*),
     &  ntmat_,ipkon(*),nactdof(0:3,*),nodeprint(*),nodeboun(*),
     &  ndirboun(*),ipompc(*),nodempc(3,*),nelem,ikboun(*),ilboun(*),
     &  ncocon(2,*)
!
      integer nk,ne,neprint,mattyp,ithermal,iprestr,i,j,k,m1,m2,jj,
     &  i1,m3,m4,kk,iener,indexe,nope,norien,iperturb,iout,noprint,
     &  nal,icmd,ihyper,nboun,nmpc,nmethod,ist,ndir,node,index,
     &  neq,kode,imat,mint3d,mint_,nopev,nfield,ndim,iorien,ielas,
     &  istiff,icmdl,ncmat_,nstate_,incrementalmpc,lb,l,jmin,jmax
!
      integer nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_
!
      real*8 co(3,*),v(0:3,*),shp(4,20),stiini(6,mint_,*),
     &  stx(6,mint_,*),stn(6,*),xl(3,20),vl(0:3,20),stre(9),
     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),
     &  alcon(0:6,ntmat_,*),vini(0:3,*),
     &  alzero(*),orab(7,*),elas(21),alph(6),rho,f(*),fn(0:3,*),
     &  skl(3,3),anisox(3,3,3,3),beta(6),q(0:3,20),vkl(0:3,3),
     &  t0(*),t1(*),prestr(6,*),eei(6,mint_,*),een(6,*),ckl(3,3),
     &  vold(0:3,*),b(*),xboun(*),coefmpc(*),vmax,eloc(9),veold(0:3,*),
     &  accold(0:3,*),elconloc(21),eth(6),xkl(3,3),voldl(0:3,20),epn(*),
     &  xikl(3,3),ener(mint_,*),enern(*),sti(6,mint_,*),
     &  eeiini(6,mint_,*),enerini(mint_,*),cocon(0:6,ntmat_,*)
!
      real*8 e,un,al,um,am1,am2,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
     &  xsj,qa,vj,t0l,t1l,bet,gam,dtime,forcempc,scal1,scal2,bnac,
     &  fixed_disp,weight,pgauss(3),vij,xkapp(6)
!
      real*8 plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
     &  xstiff(21,mint_,*),xstate(nstate_,mint_,*),plconloc(82),
     &  vokl(3,3),xstateini(nstate_,mint_,*),vikl(3,3),
     &  xstaten(nstate_,*)
!
      real*8 gauss2d1(2,1),gauss2d2(2,4),gauss2d3(2,9),gauss2d4(2,1),
     &  gauss2d5(2,3),gauss3d1(3,1),gauss3d2(3,8),gauss3d3(3,27),
     &  gauss3d4(3,1),gauss3d5(3,4),gauss3d6(3,15),gauss3d7(3,2),
     &  gauss3d8(3,9),gauss3d9(3,18),weight2d1(1),weight2d2(4),
     &  weight2d3(9),weight2d4(1),weight2d5(3),weight3d1(1),
     &  weight3d2(8),weight3d3(27),weight3d4(1),weight3d5(4),
     &  weight3d6(15),weight3d7(2),weight3d8(9),weight3d9(18)
!
      include "gauss.f"
!
      if(ithermal.le.1) then
         jmin=1
         jmax=3
      elseif(ithermal.eq.2) then
         jmin=0
         jmax=0
      else
         jmin=0
         jmax=3
      endif
!
      if((iout.ne.2).and.(iout.ne.-1)) then
!
         if((nmethod.ne.4).or.(iperturb.le.1)) then
!
!     extracting the displacement information from the solution
!
            do i=1,nk
               do j=jmin,jmax
                  if(nactdof(j,i).ne.0) then
                     v(j,i)=b(nactdof(j,i))
                  else
                     v(j,i)=0.d0
                  endif
               enddo
            enddo
!
!     for static perturbation steps v represents the incremental
!     displacements. For the total displacement vold must be added.
!
            if((iperturb.ne.0).and.(nmethod.eq.1)) then
               vmax=0.d0
               do i=1,nk
                  do j=jmin,jmax
                     vmax=max(vmax,dabs(v(j,i)))
c                     write(*,*) i,j,vold(j,i),v(j,i),vmax
                     v(j,i)=v(j,i)+vold(j,i)
                  enddo
               enddo
            endif
!
         else
!
!           direct integration dynamic step
!           b contains the acceleration increment
!
            if((ithermal.le.1).or.(ithermal.eq.3)) then
               scal1=bet*dtime*dtime
               scal2=gam*dtime
               vmax=0.d0
               do i=1,nk
                  do j=1,3
                     if(nactdof(j,i).ne.0) then
                        bnac=b(nactdof(j,i))
                        v(j,i)=vold(j,i)+scal1*bnac
                        vmax=max(vmax,dabs(scal1*bnac))
                        veold(j,i)=veold(j,i)+scal2*bnac
                        accold(j,i)=accold(j,i)+bnac
                     endif
                  enddo
               enddo
            endif
            if(ithermal.ge.2) then
               do i=1,nk
                  if(nactdof(0,i).ne.0) then
                     bnac=b(nactdof(0,i))
                     v(0,i)=vold(0,i)+bnac
                     vmax=max(vmax,dabs(bnac))
                     veold(0,i)=veold(0,i)+bnac/dtime
                  endif
               enddo
            endif
         endif
!
      endif
!
!        initialization
!
      calcul_fn=.false.
      calcul_f=.false.
      calcul_qa=.false.
      calcul_cauchy=.false.
!

!     determining which quantities have to be calculated
!
      if((iperturb.ge.2).or.((iperturb.eq.0).and.(iout.lt.0))) then
         if(iout.lt.1) then
            calcul_fn=.true.
            calcul_f=.true.
            calcul_qa=.true.
         else    
            calcul_cauchy=.true.
         endif
      endif
      if(((noelplab(5).eq.'RF').or.(nodeflab(5).eq.'RF'))
     &     .and.(iout.gt.0)) calcul_fn=.true.
!
!     initializing fn
!
      if(calcul_fn) then
         do i=1,nk
            do j=1,3
               fn(j,i)=0.d0
            enddo
         enddo
      endif
!
!     initializing f
!
      if(calcul_f) then
         do i=1,neq
            f(i)=0.d0
         enddo
      endif
!
!     inserting the boundary conditions
!
      do i=1,nboun
         if(ndirboun(i).gt.3) cycle
         fixed_disp=xboun(i)
         if((nmethod.eq.4).and.(iperturb.gt.1)) then
            ndir=ndirboun(i)
            node=nodeboun(i)
            accold(ndir,node)=(xboun(i)-v(ndir,node))/
     &        (bet*dtime*dtime)
            veold(ndir,node)=veold(ndir,node)+
     &        gam*dtime*accold(ndir,node)
         endif
         v(ndirboun(i),nodeboun(i))=fixed_disp
c         write(*,*) 'boun ',nodeboun(i),ndirboun(i),fixed_disp
      enddo
!
!     inserting the mpc information
!     the parameter incrementalmpc indicates whether the
!     incremental displacements enter the mpc or the total 
!     displacements (incrementalmpc=0)
!
      do i=1,nmpc
         if((labmpc(i)(1:20).eq.'                    ').or.
     &      (labmpc(i)(1:6).eq.'CYCLIC').or.
     &      (labmpc(i)(1:5).eq.'RIGID').or.
     &      (labmpc(i)(1:9).eq.'SUBCYCLIC')) then
            incrementalmpc=0
         else
            incrementalmpc=1
         endif
         ist=ipompc(i)
         node=nodempc(1,ist)
         ndir=nodempc(2,ist)
         index=nodempc(3,ist)
         fixed_disp=0.d0
         if(index.ne.0) then
            do
               if(incrementalmpc.eq.0) then
                  fixed_disp=fixed_disp-coefmpc(index)*
     &                 v(nodempc(2,index),nodempc(1,index))
               else
c                  write(*,102) 
c     &             nodempc(1,index),nodempc(2,index),
c     &             v(nodempc(2,index),nodempc(1,index)),
c     &             vold(nodempc(2,index),nodempc(1,index))
c 102              format('results,node,dir,newvalue,oldvalue',2i5,
c     &                2(1x,e11.4))
                  fixed_disp=fixed_disp-coefmpc(index)*
     &                 (v(nodempc(2,index),nodempc(1,index))-
     &                  vold(nodempc(2,index),nodempc(1,index)))
c                  write(*,*) fixed_disp,
c     & vold(nodempc(2,index),nodempc(1,index)),
c     & v(nodempc(2,index),nodempc(1,index))
               endif
               index=nodempc(3,index)
               if(index.eq.0) exit
            enddo
         endif
         fixed_disp=fixed_disp/coefmpc(ist)
c                  write(*,*) fixed_disp
         if(incrementalmpc.eq.1) then
            fixed_disp=fixed_disp+vold(ndir,node)
c            write(*,*) 'results ',i,vold(ndir,node),fixed_disp
         endif
         if((nmethod.eq.4).and.(iperturb.gt.1)) then
            accold(ndir,node)=(fixed_disp-v(ndir,node))/
     &        (bet*dtime*dtime)
            veold(ndir,node)=veold(ndir,node)+
     &        gam*dtime*accold(ndir,node)
         endif
         v(ndir,node)=fixed_disp
c         write(*,102) node,ndir,v(ndir,node),vold(ndir,node)
      enddo
!
!     check whether there are any strain output requests
!
      if ((noelplab(7).eq.'ENER').or.(nodeflab(7).eq.'ENER')) then
         iener=1
      else
         iener=0
      endif
!
      qa=0.d0
      nal=0
!
!     calculation of the stresses in the integration points
! 
      if((ithermal.le.1).or.(ithermal.eq.3)) then
!
      do i=1,ne
!
         if(ipkon(i).lt.0) cycle
         imat=ielmat(i)
         amat=matname(imat)
         if(norien.gt.0) then
            iorien=ielorien(i)
         else
            iorien=0
         endif
!
c         if(ipkon(i).lt.0) cycle
         indexe=ipkon(i)
         if(lakon(i)(4:4).eq.'2') then
            nope=20
         elseif(lakon(i)(4:4).eq.'8') then
            nope=8
         elseif(lakon(i)(4:5).eq.'10') then
            nope=10
         elseif(lakon(i)(4:4).eq.'4') then
            nope=4
         elseif(lakon(i)(4:5).eq.'15') then
            nope=15
         else
            nope=6
         endif
!
         if(lakon(i)(4:5).eq.'8R') then
            mint3d=1
         elseif((lakon(i)(4:4).eq.'8').or.
     &          (lakon(i)(4:6).eq.'20R')) then
            mint3d=8
         elseif(lakon(i)(4:4).eq.'2') then
            mint3d=27
         elseif(lakon(i)(4:5).eq.'10') then
            mint3d=4
         elseif(lakon(i)(4:4).eq.'4') then
            mint3d=1
         elseif(lakon(i)(4:5).eq.'15') then
            mint3d=9
         else
            mint3d=2
         endif
!
         do j=1,nope
            konl(j)=kon(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl(j))
               vl(k,j)=v(k,konl(j))
               voldl(k,j)=vold(k,konl(j))
            enddo
         enddo
!
!        check for hyperelastic material
!
         if(nelcon(1,imat).lt.0) then
            ihyper=1
         else
            ihyper=0
         endif
!
!        q contains the nodal forces per element; initialisation of q
!
         if((iperturb.ge.2).or.((iperturb.eq.0).and.(iout.lt.1))) then
            do m1=1,nope
               do m2=1,3
                  q(m2,m1)=fn(m2,konl(m1))
               enddo
            enddo
         endif
!
         do jj=1,mint3d
            if(lakon(i)(4:5).eq.'8R') then
               xi=gauss3d1(1,jj)
               et=gauss3d1(2,jj)
               ze=gauss3d1(3,jj)
               weight=weight3d1(jj)
            elseif((lakon(i)(4:4).eq.'8').or.
     &             (lakon(i)(4:6).eq.'20R'))
     &        then
               xi=gauss3d2(1,jj)
c               if(nope.eq.20) xi=xi+1.d0
               et=gauss3d2(2,jj)
               ze=gauss3d2(3,jj)
               weight=weight3d2(jj)
            elseif(lakon(i)(4:4).eq.'2') then
c               xi=gauss3d3(1,jj)+1.d0
               xi=gauss3d3(1,jj)
               et=gauss3d3(2,jj)
               ze=gauss3d3(3,jj)
               weight=weight3d3(jj)
            elseif(lakon(i)(4:5).eq.'10') then
               xi=gauss3d5(1,jj)
               et=gauss3d5(2,jj)
               ze=gauss3d5(3,jj)
               weight=weight3d5(jj)
            elseif(lakon(i)(4:4).eq.'4') then
               xi=gauss3d4(1,jj)
               et=gauss3d4(2,jj)
               ze=gauss3d4(3,jj)
               weight=weight3d4(jj)
            elseif(lakon(i)(4:5).eq.'15') then
               xi=gauss3d8(1,jj)
               et=gauss3d8(2,jj)
               ze=gauss3d8(3,jj)
               weight=weight3d8(jj)
            else
               xi=gauss3d7(1,jj)
               et=gauss3d7(2,jj)
               ze=gauss3d7(3,jj)
               weight=weight3d7(jj)
            endif
!
            if(nope.eq.20) then
               call shape20h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp)
            else
               call shape6w(xi,et,ze,xl,xsj,shp)
            endif
!
!                 vkl(m2,m3) contains the derivative of the m2-
!                 component of the displacement with respect to
!                 direction m3
!
            do m2=1,3
               do m3=1,3
                  vkl(m2,m3)=0.d0
               enddo
            enddo
!
            do m1=1,nope
               do m2=1,3
                  do m3=1,3
                     vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
                  enddo
               enddo
            enddo
!
            kode=nelcon(1,imat)
!
!                 calculating the strain
!
!                 attention! exy,exz and eyz are engineering strains!
!
            exx=vkl(1,1)
            eyy=vkl(2,2)
            ezz=vkl(3,3)
            exy=vkl(1,2)+vkl(2,1)
            exz=vkl(1,3)+vkl(3,1)
            eyz=vkl(2,3)+vkl(3,2)
!
!                 Lagrangian/Eulerian strain tensor for geometric
!                 nonlinear calculations
!
            if(iperturb.ge.2) then
!
!                    Lagrangian strain
!
               exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
               eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
               ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
               exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
     &              vkl(3,1)*vkl(3,2)
               exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
     &              vkl(3,1)*vkl(3,3)
               eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
     &              vkl(3,2)*vkl(3,3)
!
            endif
!
!                 storing the local strains
!
            if((iperturb.ne.0).or.(iout.ge.0)) then
               eloc(1)=exx
               eloc(2)=eyy
               eloc(3)=ezz
               eloc(4)=exy/2.d0
               eloc(5)=exz/2.d0
               eloc(6)=eyz/2.d0
            else
!
!                    linear iteration within a nonlinear increment:
!                    tangent matrix is calculated at strain at
!                    the end of the previous increment
!
               do m1=1,6
                  eloc(m1)=eei(m1,jj,i)
               enddo
            endif
!
!                 calculating the deformation gradient (needed to
!                 convert the element stiffness matrix from spatial
!                 coordinates to material coordinates
!                 deformation plasticity)
!
            if((kode.eq.-50).or.(kode.le.-100)) then
!
!                    calculating the deformation gradient
!
               xkl(1,1)=vkl(1,1)+1
               xkl(2,2)=vkl(2,2)+1.
               xkl(3,3)=vkl(3,3)+1.
               xkl(1,2)=vkl(1,2)
               xkl(1,3)=vkl(1,3)
               xkl(2,3)=vkl(2,3)
               xkl(2,1)=vkl(2,1)
               xkl(3,1)=vkl(3,1)
               xkl(3,2)=vkl(3,2)
!
!                    calculating the Jacobian
!
               vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
     &              -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
     &              +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
!
!              inversion of the deformation gradient (only for
!              deformation plasticity)
!
               if(kode.eq.-50) then
!
                  ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
                  ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
                  ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
                  ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
                  ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
                  ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
                  ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
                  ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
                  ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
!
!                 converting the Lagrangian strain into Eulerian
!                 strain (only for deformation plasticity)
!
                  cauchy=.false.
                  call str2mat(eloc,ckl,vj,cauchy)
               endif
!                        
            endif
!
!                 calculating fields for incremental plasticity
!
            if(kode.le.-100) then
!
!              calculating the deformation gradient at the
!              start of the increment
!
!              calculating the displacement gradient at the
!              start of the increment
!
               do m2=1,3
                  do m3=1,3
                     vikl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vikl(m2,m3)=vikl(m2,m3)
     &                       +shp(m3,m1)*vini(m2,konl(m1))
                     enddo
                  enddo
               enddo
!
!              calculating the deformation gradient of the old
!              fields
!
               xikl(1,1)=vikl(1,1)+1
               xikl(2,2)=vikl(2,2)+1.
               xikl(3,3)=vikl(3,3)+1.
               xikl(1,2)=vikl(1,2)
               xikl(1,3)=vikl(1,3)
               xikl(2,3)=vikl(2,3)
               xikl(2,1)=vikl(2,1)
               xikl(3,1)=vikl(3,1)
               xikl(3,2)=vikl(3,2)
!
!              calculating the Jacobian
!
               vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,2))
     &              -xikl(1,2)*(xikl(2,1)*xikl(3,3)
     &              -xikl(2,3)*xikl(3,1))
     &              +xikl(1,3)*(xikl(2,1)*xikl(3,2)
     &              -xikl(2,2)*xikl(3,1))
!
!              stresses at the start of the increment
!
               do m1=1,6
                  stre(m1)=stiini(m1,jj,i)
               enddo
!
            endif
!
!                 prestress values
!
            if(iprestr.eq.0) then
               do kk=1,6
                  beta(kk)=0.d0
               enddo
            else
               do kk=1,6
                  beta(kk)=-prestr(kk,i)
               enddo
            endif
!
            if(ithermal.ge.1) then
!
!              calculating the temperature difference in
!              the integration point
!
               t0l=0.d0
               t1l=0.d0
               if(ithermal.eq.1) then
                  if(lakon(i)(4:5).eq.'8 ') then
                     do i1=1,nope
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+t1(konl(i1))/8.d0
                     enddo
                  elseif(lakon(i)(4:6).eq.'20 ') then
                     call lintemp(t0,t1,konl,nope,jj,t0l,t1l)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*t1(konl(i1))
                     enddo
                  endif
               elseif(ithermal.ge.2) then
                  if(lakon(i)(4:5).eq.'8 ') then
                     do i1=1,nope
                        t0l=t0l+t0(konl(i1))/8.d0
                        t1l=t1l+vold(0,konl(i1))/8.d0
                     enddo
                  elseif(lakon(i)(4:6).eq.'20 ') then
                     call lintemp_th(t0,vold,konl,nope,jj,t0l,t1l)
                  else
                     do i1=1,nope
                        t0l=t0l+shp(4,i1)*t0(konl(i1))
                        t1l=t1l+shp(4,i1)*vold(0,konl(i1))
                     enddo
                  endif
               endif
               tt=t1l-t0l
            endif
!
!                 calculating the coordinates of the integration point
!                 for material orientation purposes (for cylindrical
!                 coordinate systems)
!
            if(iorien.gt.0) then
               do j=1,3
                  pgauss(j)=0.d0
                  do i1=1,nope
                     pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
                  enddo
               enddo
            endif
!
!                 material data; for linear elastic materials
!                 this includes the calculation of the stiffness
!                 matrix
!
            istiff=0
!
            call materialdata(elcon,nelcon,rhcon,nrhcon,alcon,
     &           nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
     &           elas,alph,rho,i,ithermal,alzero,mattyp,
     &           t0l,t1l,ihyper,istiff,
     &           elconloc,eth,kode,plicon,nplicon,
     &           plkcon,nplkcon,npmat_,
     &           plconloc,mint_,dtime,i,jj,
     &           xstiff,ncmat_)
!
            if((ihyper.eq.0).and.
     &           ((iperturb.ne.0).or.(iout.ge.0))) then
               if(mattyp.eq.1) then
                  e=elas(1)
                  un=elas(2)
                  al=un*e/(1.d0+un)/(1.d0-2.d0*un)
                  um=e/2.d0/(1.d0+un)
                  am1=al+2.d0*um
                  am2=2.d0*um
               elseif(mattyp.eq.3) then
                  call anisotropic(elas,anisox)
               endif
!
               icmdl=3
               call linel(ithermal,mattyp,beta,al,um,am1,alph,tt,
     &              elas,icmdl,exx,eyy,ezz,exy,exz,eyz,stre,
     &              anisox)
!
            elseif(ihyper.ne.0) then
               mattyp=3
               call hyperel(elconloc,elas,eloc,kode,eth,ithermal,
     &              icmd,beta,stre,xkl,ckl,vj,xikl,vij,
     &              plconloc,xstate,xstateini,ielas,
     &              amat,t1l,dtime,i,jj,nstate_,mint_,
     &              iorien,pgauss,orab)
               do m1=1,21
                  xstiff(m1,jj,i)=elas(m1)
               enddo
            endif
!
            if((iperturb.eq.0).and.(iout.lt.0)) then
!
!                    if the forced displacements were changed at
!                    the start of a nonlinear step, the nodal
!                    forces due do this displacements are 
!                    calculated in a purely linear way, and
!                    the first iteration is purely linear in order
!                    to allow the displacements to redistribute
!                    in a quasi-static way (only applies to
!                    quasi-static analyses (*STATIC))
!
               do m2=1,3
                  do m3=1,3
                     vokl(m2,m3)=0.d0
                  enddo
               enddo
!
               do m1=1,nope
                  do m2=1,3
                     do m3=1,3
                        vokl(m2,m3)=vokl(m2,m3)+
     &                       shp(m3,m1)*voldl(m2,m1)
                     enddo
                  enddo
               enddo
!
               eloc(1)=exx-vokl(1,1)
               eloc(2)=eyy-vokl(2,2)
               eloc(3)=ezz-vokl(3,3)
               eloc(4)=exy-(vokl(1,2)+vokl(2,1))
               eloc(5)=exz-(vokl(1,3)+vokl(3,1))
               eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
               if(mattyp.eq.1) then
                  e=elas(1)
                  un=elas(2)
                  um=e/(1.d0+un)
                  al=un*um/(1.d0-2.d0*un)
                  um=um/2.d0
                  am1=al*(eloc(1)+eloc(2)+eloc(3))
                  stre(1)=am1+2.d0*um*eloc(1)
                  stre(2)=am1+2.d0*um*eloc(2)
                  stre(3)=am1+2.d0*um*eloc(3)
                  stre(4)=um*eloc(4)
                  stre(5)=um*eloc(5)
                  stre(6)=um*eloc(6)
               elseif(mattyp.eq.2) then
                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
     &                 +eloc(3)*elas(4)
                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
     &                 +eloc(3)*elas(5)
                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
     &                 +eloc(3)*elas(6)
                  stre(4)=eloc(4)*elas(7)
                  stre(5)=eloc(5)*elas(8)
                  stre(6)=eloc(6)*elas(9)
               elseif(mattyp.eq.3) t

454-calculix

Exec Time : 1.306602e+08 cycles

Matching : 40.0%

__invivo__SPOOLES_Chv_src_Chv_factor_symmetric1x1_600

Real time: 1.753996e+06 -- Predicted time: 2.873963e+06

Callcount: 720

__invivo__mastruct_mastruct_615

Real time: 4.825128e+06 -- Predicted time: 4.780256e+06

Callcount: 1

__invivo__mastruct_mastruct_331

Real time: 1.939523e+07 -- Predicted time: 1.919306e+07

Callcount: 1

__invivo__SPOOLES_Utilities_src_Utilities_DV_DVdot33_1241

Real time: 3.388256e+06 -- Predicted time: 7.024252e+06

Callcount: 7920

__invivo__mastruct_mastruct_705

Real time: 4.270208e+06 -- Predicted time: 4.277338e+06

Callcount: 1

__invivo__SPOOLES_Utilities_src_Utilities_newsort_IVqsortUp_975

Real time: 1.825840e+06 -- Predicted time: 1.007780e+07

Callcount: 12588

__invivo__isortii_isortii__199

Real time: 1.202479e+07 -- Predicted time: 9.773259e+06

Callcount: 1441

__invivo__results_results__1301

Real time: 4.682592e+06 -- Predicted time: 1.284400e+05

Callcount: 1

__invivo__SPOOLES_Utilities_src_Utilities_sortAndCompress_IV2DVsortUpAndCompress_339

Real time: 3.298192e+06 -- Predicted time: 3.343572e+06

Callcount: 1

invivoSPOOLES_Chv_src_Chv_factor_symmetric1x1_600

invivomastruct_mastruct_615

invivomastruct_mastruct_331

invivoSPOOLES_Utilities_src_Utilities_DV_DVdot33_1241

invivomastruct_mastruct_705

invivoSPOOLES_Utilities_src_Utilities_newsort_IVqsortUp_975

invivoisortii_isortii__199

invivoresults_results__1301

invivoSPOOLES_Utilities_src_Utilities_sortAndCompress_IV2DVsortUpAndCompress_339