void *
btreeCreateNodeBlock(GdbBlock *block, void *extra)
{
    BTreeNode *node;
    BTree *tree;

    tree = (BTree *)extra;

    MEM_CHECK(node = (BTreeNode *)SAFE_MALLOC(sizeof(BTreeNode), LOC_BTREE_0166));
    memset(node, 0, sizeof(BTreeNode));

    node->tree  = tree;
    node->block = block;

    /*comment: when reach here, we have no idea about keyCount, hence alloc children/keySizes/keys as the max possible num: the order*/

    MEM_CHECK(node->children = (offset_t *)SAFE_MALLOC(tree->order * sizeof(offset_t), LOC_BTREE_0167));
    memset(node->children, 0, tree->order * sizeof(offset_t));


    MEM_CHECK(node->keySizes = (uint16_t *)SAFE_MALLOC((tree->order - 1) * sizeof(uint16_t), LOC_BTREE_0168));
    memset(node->keySizes, 0, (tree->order - 1)  * sizeof(uint16_t));


    MEM_CHECK(node->keys = (uint8_t **)SAFE_MALLOC((tree->order - 1) * sizeof(uint8_t *), LOC_BTREE_0169));
    memset(node->keys, 0, (tree->order - 1) * sizeof(uint8_t *));

    return node;
}

void PetscMatrix::alloc() {
  _F_
#ifdef WITH_PETSC
  assert(pages != NULL);

  // calc nnz
  int *nnz_array = new int[size];
  MEM_CHECK(nnz_array);

  // fill in nnz_array
  int aisize = get_num_indices();
  int *ai = new int[aisize];
  MEM_CHECK(ai);

  // sort the indices and remove duplicities, insert into ai
  int pos = 0;
  for (unsigned int i = 0; i < size; i++) {
    nnz_array[i] = sort_and_store_indices(pages[i], ai + pos, ai + aisize);
    pos += nnz_array[i];
  }
  // stote the number of nonzeros
  nnz = pos;
  delete [] pages; pages = NULL;
  delete [] ai;

  //
  MatCreateSeqAIJ(PETSC_COMM_SELF, size, size, 0, nnz_array, &matrix);
//	MatSetOption(matrix, MAT_ROW_ORIENTED);
//	MatSetOption(matrix, MAT_ROWS_SORTED);

  delete [] nnz_array;

  inited = true;
#endif
}

void SuperLUMatrix::alloc()
{
  _F_
  assert(pages != NULL);
  assert(size > 0);
  
  // Initialize the arrays Ap and Ai.
  Ap = new int [size + 1];
  MEM_CHECK(Ap);
  int aisize = get_num_indices();
  Ai = new int [aisize];
  MEM_CHECK(Ai);
  
  // sort the indices and remove duplicities, insert into Ai
  int i, pos = 0;
  for (i = 0; i < size; i++) {
    Ap[i] = pos;
    pos += sort_and_store_indices(pages[i], Ai + pos, Ai + aisize);
  }
  Ap[i] = pos;
  
  delete [] pages;
  pages = NULL;
  
  nnz = Ap[size];

  Ax = new slu_scalar [nnz];
  memset(Ax, 0, sizeof(slu_scalar) * nnz);
}

void UMFPackMatrix::alloc() {
	_F_
	assert(pages != NULL);

	// initialize the arrays Ap and Ai
	Ap = new int [size + 1];
	MEM_CHECK(Ap);
	int aisize = get_num_indices();
	Ai = new int [aisize];
	MEM_CHECK(Ai);

	// sort the indices and remove duplicities, insert into Ai
	int i, pos = 0;
	for (i = 0; i < size; i++) {
		Ap[i] = pos;
		pos += sort_and_store_indices(pages[i], Ai + pos, Ai + aisize);
	}
	Ap[i] = pos;

	delete [] pages;
	pages = NULL;

	Ax = new scalar [Ap[size]];
	MEM_CHECK(Ax);
	memset(Ax, 0, sizeof(scalar) * Ap[size]);
}

void CSCMatrix::alloc() {
  _F_
  assert(pages != NULL);
  if (size <= 0)
      error("UMFPack failed, matrix size must be greater than 0.");

  // initialize the arrays Ap and Ai
  Ap = new int [size + 1];
  MEM_CHECK(Ap);
  int aisize = get_num_indices();
  Ai = new int [aisize];
  MEM_CHECK(Ai);

  // sort the indices and remove duplicities, insert into Ai
  unsigned int i;
  int pos = 0;
  for (i = 0; i < size; i++) {
    Ap[i] = pos;
    pos += sort_and_store_indices(pages[i], Ai + pos, Ai + aisize);
  }
  Ap[i] = pos;

  delete [] pages;
  pages = NULL;

  nnz = Ap[size];
  
  Ax = new scalar [nnz];
  MEM_CHECK(Ax);
  memset(Ax, 0, sizeof(scalar) * nnz);
}

void RefMap::calc_face_normal(int iface, const int np, const QuadPt3D *pt, double *&nx, double *&ny, double *&nz) {
	_F_
	assert(mesh != NULL);

	double3x3 *m = get_ref_map(np, pt);

	nx = new double[np]; MEM_CHECK(nx);
	ny = new double[np]; MEM_CHECK(ny);
	nz = new double[np]; MEM_CHECK(nz);

	// FIXME: refactor this!
	// - it would be nice if calculation of normals would be based on the same algorithm for all elements
	//   e.g. to take a normal from ref. domain and transform it to the physical one
	int t_dir_1, t_dir_2; //directions of tangents ot the reference face such that t_dir_1 x t_dir_2 = outer normal
	switch (element->get_mode()) {
		case MODE_TETRAHEDRON: {
			const int *face_vtx = element->get_face_vertices(iface);
			Vertex vtx[Tri::NUM_VERTICES];
			for (int i = 0; i < element->get_num_vertices(); i++)
				vtx[i] = vertex[face_vtx[i]];

			Point3D v1 = { vtx[1].x - vtx[0].x, vtx[1].y - vtx[0].y, vtx[1].z - vtx[0].z };
			Point3D v2 = { vtx[2].x - vtx[0].x, vtx[2].y - vtx[2].y, vtx[2].z - vtx[0].z };
			Point3D n = normalize(cross_product(v1, v2));

			for (int i = 0; i < np; i++) {
				nx[i] = n.x;
				ny[i] = n.y;
				nz[i] = n.z;
			}
			} break;

		case MODE_HEXAHEDRON:
			switch (iface) {
				case 0: t_dir_1 = 2; t_dir_2 = 1; break;
				case 1: t_dir_1 = 1; t_dir_2 = 2; break;
				case 2: t_dir_1 = 0; t_dir_2 = 2; break;
				case 3: t_dir_1 = 2; t_dir_2 = 0; break;
				case 4: t_dir_1 = 1; t_dir_2 = 0; break;
				case 5: t_dir_1 = 0; t_dir_2 = 1; break;
			}
			for (int i = 0; i < np; i++) {
				Point3D tangent1 = { m[i][0][t_dir_1], m[i][1][t_dir_1], m[i][2][t_dir_1] };
				Point3D tangent2 = { m[i][0][t_dir_2], m[i][1][t_dir_2], m[i][2][t_dir_2] };
				Point3D normal = normalize(cross_product(tangent1, tangent2));

				nx[i] = normal.x;
				ny[i] = normal.y;
				nz[i] = normal.z;
			}
			break;

		case MODE_PRISM:
			EXIT(HERMES_ERR_NOT_IMPLEMENTED);
	}

	delete [] m;
}

void EpetraMatrix::prealloc(unsigned int n)
{
  _F_
#ifdef HAVE_EPETRA
  this->size = n;
  // alloc trilinos structs
  std_map = new Epetra_Map(n, 0, seq_comm); MEM_CHECK(std_map);
  grph = new Epetra_CrsGraph(Copy, *std_map, 0); MEM_CHECK(grph);
#endif
}

void EpetraMatrix::alloc()
{
  _F_
#ifdef HAVE_EPETRA
  grph->FillComplete();
  // create the matrix
  mat = new Epetra_CrsMatrix(Copy, *grph); MEM_CHECK(mat);
#ifdef HERMES_COMMON_COMPLEX
  mat_im = new Epetra_CrsMatrix(Copy, *grph); MEM_CHECK(mat_im);
#endif
#endif
}

/************************************************************************* 
 * f_sql_handles( void ) 
 *
 * returns an array containing the IDs of the currently opened connections
 *************************************************************************/
svalue_t *
f_sql_handles( svalue_t * argv, int argc )
{
#if ODBC_DEBUG & DEBUG_FUNC
   printf( "call f_sql_handles( )\n" );
#endif
   int  cnt, i;
   hDBC * tmp;

   vector_t    * vec;

   argv++;
   
   if ( !hODBCEnv )
      init_odbc_environment();
  
   if ( !hODBCEnv ) {
      put_number( argv, 0 );
      return( argv );
   }
   
   if ( !hODBCEnv->hDBCons ) {

      vec = allocate_array( 0 );
      MEM_CHECK( vec );

      put_array( argv, vec );
      return( argv );
   }

   tmp = hODBCEnv->hDBCons;
   cnt = 1;

   while( (tmp = tmp->next ) )
      cnt++;

   vec = allocate_array( cnt );
   MEM_CHECK( vec );

   tmp = hODBCEnv->hDBCons;
   for ( i = 0; i < cnt; ++i ) {
      put_number( vec->item + i, tmp->ID );
      tmp = tmp->next;
   }

   put_array( argv, vec );
#if ODBC_DEBUG & DEBUG_FUNC
   printf( "ret f_sql_handles( )\n" );
#endif
   return( argv );
}

static mapping_t * 
fetch_into_mapping( hDBC * handle )
{
   int       i;

   mapping_t  * map;   
   svalue_t  * key,
             * value;

   STORE_DOUBLE_USED;

   map = allocate_mapping( handle->colcnt, 1 );
   MEM_CHECK( map );
         
   for( i = 0; i < handle->colcnt; ++i ) {
      if ( !handle->columns[ i ] ) continue;

      //printf( " fetch_into_mapping[%2d] ", i );
      
      key = pxalloc( sizeof( svalue_t ) ); 
      MEM_CHECK( key );
         
      put_malloced_string( key, string_copy( handle->columns[ i ]->name ) );
         
      value = get_map_lvalue( map, key );
      MEM_CHECK( value );

      switch( handle->columns[ i ]->type ){
         case T_FLOAT:
            //printf( "float=%f\n",  handle->columns[ i ]->data.double_v );
            STORE_DOUBLE( value, *handle->columns[ i ]->data.double_v );

            value->type = T_FLOAT;
            break; 

         case T_NUMBER:
            //printf( "number=%d\n",  *handle->columns[ i ]->data.number_v );
            put_number( value, *handle->columns[ i ]->data.number_v );
            break;

         case T_STRING: 
         default      :
            //printf( "string=%s\n",  handle->columns[ i ]->data.string_v );
            put_malloced_string( value, string_copy(  handle->columns[ i ]->data.string_v ) );
            break;
      }
   }

   return( map );
}

RawFile * rawFileOpen(const uint8_t *file_name, const int flags, const uint32_t file_size, const word_t cdfs_md_id)
{
    MEM_CHECK(file_name);

    if(flags & O_RDWR)
    {
        RawFile   * raw_file;
#if 0
        CBYTES    * cbytes;
        CSTRING   * fname_cstr;
#endif
        raw_file = rawFileNew(file_name, -1, O_RDWR, file_size, cdfs_md_id);
        if(NULL == raw_file)
        {
            dbg_log(SEC_0132_RAW, 0)(LOGSTDOUT, "error:rawFileOpen: new raw file failed\n");
            return NULL;
        }

        if(RAW_FILE_SUCC != rawFileLoad(raw_file))
        {
            dbg_log(SEC_0132_RAW, 0)(LOGSTDOUT, "error:rawFileOpen: load %s failed\n", (char *)file_name);
            raw_file->fd = -1;
            rawFileFree(raw_file);
            return NULL;
        }
        return raw_file;
    }

    if(flags & O_CREAT)
    {
        return rawFileCreate(file_name, file_size, cdfs_md_id);
    }

    return NULL;
}

double *RefMap::get_jacobian(const int np, const QuadPt3D *pt, bool trans) {
	_F_

	double *jac = new double[np]; MEM_CHECK(jac);
	if (is_const_jacobian) {
		if (trans)
			for (int i = 0; i < np; i++)
				jac[i] = const_jacobian * pt[i].w;
		else
			for (int i = 0; i < np; i++)
				jac[i] = const_jacobian;
	}
	else {
		double3x3 *m = get_ref_map(np, pt);

		double trj = get_transform_jacobian();
		if (trans)
			for (int i = 0; i < np; i++)
				jac[i] = det(m[i]) * trj * pt[i].w;
		else
			for (int i = 0; i < np; i++)
				jac[i] = det(m[i]) * trj;

		delete [] m;
	}

	return jac;
}

static uint8_t *__dupFileName(const uint8_t *root_path)
{
    uint8_t *file_name;
    MEM_CHECK(file_name = (uint8_t *)SAFE_MALLOC(strlen((char *)root_path) + 1, LOC_RAW_0001));
    sprintf((char *)file_name, "%s", (char *)root_path);
    return (file_name);
}

uint8_t   rawFileInit(RawFile *raw_file, const uint8_t *file_name, const int fd, const int flags, const uint32_t file_size, const word_t cdfs_md_id)
{
    RawData *raw_data;
    MEM_CHECK(raw_file);
    raw_data = rawDataNew(file_size, RAW_FILE_HEAD_SIZE);
    if(NULL == raw_data)
    {
        dbg_log(SEC_0132_RAW, 0)(LOGSTDOUT,"error:rawFileInit: new raw data failed\n");
        return RAW_FILE_FAIL;
    }
    raw_file->fd = fd;
    raw_file->open_flags = flags;
    raw_file->cdfs_md_id = cdfs_md_id;

    if(NULL != file_name)
    {
        raw_file->file_name = __dupFileName(file_name);
    }
    else
    {
        raw_file->file_name = NULL;
    }
    raw_file->raw_data  = raw_data;

    return RAW_FILE_SUCC;
}

uint8_t   rawFileUpdate8s(RawFile *raw_file, const uint8_t *data, const uint32_t len, const uint32_t offset)
{
    MEM_CHECK(raw_file);
    //MEM_CHECK(offset);

    return rawDataUpdate8s(raw_file->raw_data, offset, data, len);
}

uint8_t   rawFileAppend8s(RawFile *raw_file, const uint8_t *data, const uint32_t len, uint32_t *offset)
{
    MEM_CHECK(raw_file);
    //MEM_CHECK(offset);
    (*offset) = raw_file->raw_data->cur_size;
    return rawDataPut8s(raw_file->raw_data, (*offset), data, len);
}

void
gdbWriteFreeBlockList(GDatabase *db, GdbFreeBlock *blocks, uint32_t count)
{
    uint32_t listSize;
    uint8_t *buffer;
    uint32_t i, counter = 0;
    offset_t __offset;

    if (db == NULL || blocks == NULL)
    {
        return;
    }
    /* Get the total size of the list. */
    listSize = sizeof(uint32_t) + count * (sizeof(uint16_t) + sizeof(offset_t));

    /* Allocate the buffer for the block list. */
    MEM_CHECK(buffer = (uint8_t *)SAFE_MALLOC(listSize, LOC_DB_0006));

    gdbPut32(buffer, &counter, count);

    for (i = 0; i < count; i++)
    {
        gdbPut16(buffer, &counter, blocks[i].size);
        gdbPutOffset(buffer, &counter, blocks[i].offset);
    }

    rawFileSeek(db->idxRawFile, DB_FREE_BLOCK_LIST_OFFSET, SEEK_SET);
    __offset = DB_FREE_BLOCK_LIST_OFFSET;

    rawFileWrite(db->idxRawFile, __offset, buffer, listSize, 1, LOC_DB_0007);

    SAFE_FREE(buffer, LOC_DB_0008);
}

/************************************************************************* 
 * allocate_db_connection( ) . allocates a hDBC struct with the given ID                      
 * push_db_connection( ) ..... adds an allocated hDBC struct to the DB
 *                             connections 
 * pop_db_connection( ) ...... removes an allocated hDBC struct from  
 *                             the DB connections
 * get_db_connection_by_id( )  returns the requested handle and moves it on    
 *                             top of hODBCEnv->hDBCons                               
 * dispose_db_connection( ) .. disposes the passed connection         
 *************************************************************************/
static hDBC *
allocate_db_connection( int ID )
{
   hDBC      * handle;

   if ( !ID )
      return( NULL );

   handle = pxalloc( sizeof( hDBC ) );
   MEM_CHECK( handle );
   
   handle->ID = ID;
   
   handle->name    = NULL;
   
   handle->prev    = NULL;
   handle->next    = NULL;

   handle->hDBCon  = NULL;
   handle->hStmt   = NULL;

   handle->colcnt  = 0;
   handle->rowcnt  = 0;
   handle->columns = NULL;

   return( handle );
}

GdbBlock *
gdbNewBlock(GDatabase *db, blocktype_t blockType, void *extra)
{
	GdbBlock *block;
	blocktype_t typeIndex;

	if (db == NULL || !GDB_VALID_BLOCK_TYPE(blockType))
		return NULL;

	MEM_CHECK(block = (GdbBlock *)malloc(sizeof(GdbBlock)));
	memset(block, 0, sizeof(GdbBlock));

	block->type     = blockType;
	block->db       = db;
	block->inList   = 0;
	block->refCount = 0;

	GDB_SET_DIRTY(block);

	typeIndex = blockType - 1;

	block->multiple = blockTypeInfo[typeIndex].multiple;

	if (blockTypeInfo[typeIndex].create != NULL)
	{
		block->detail = blockTypeInfo[typeIndex].create(block, extra);
	}

	return block;
}

    bool AmesosSolver<std::complex<double> >::solve()
    {
      _F_;
      assert(m != NULL);
      assert(rhs != NULL);

      assert(m->size == rhs->size);

      Hermes::TimePeriod tmr;  

      error("AmesosSolver<Scalar>::solve() not yet implemented for complex problems");

      if (!setup_factorization())
      {
        warning("AmesosSolver: LU factorization could not be completed");
        return false;
      }

      int status = solver->Solve();
      if (status != 0) 
      {
        error("AmesosSolver: Solution failed.");
        return false;
      }

      tmr.tick();
      this->time = tmr.accumulated();

      delete [] this->sln;
      this->sln = new std::complex<double>[m->size]; MEM_CHECK(this->sln);
      // copy the solution into sln vector
      memset(this->sln, 0, m->size * sizeof(std::complex<double>));

      return true;
    }