int igraph_delete_vertices(igraph_t *graph, igraph_vs_t vertices) {

  long int no_of_edges=igraph_ecount(graph);
  long int no_of_nodes=igraph_vcount(graph);
  igraph_vector_t edge_recoding, vertex_recoding;
  igraph_vit_t vit;
  igraph_t newgraph;
  long int i, j;
  long int remaining_vertices, remaining_edges;
  
  IGRAPH_VECTOR_INIT_FINALLY(&vertex_recoding, no_of_nodes);
  IGRAPH_VECTOR_INIT_FINALLY(&edge_recoding, no_of_edges);
 
  IGRAPH_CHECK(igraph_vit_create(graph, vertices, &vit));
  IGRAPH_FINALLY(igraph_vit_destroy, &vit);
  
  /* mark the vertices to delete */
  for (; !IGRAPH_VIT_END(vit); IGRAPH_VIT_NEXT(vit) ) {
    long int vertex=IGRAPH_VIT_GET(vit);
    if (vertex < 0 || vertex >= no_of_nodes) {
      IGRAPH_ERROR("Cannot delete vertices", IGRAPH_EINVVID);
    }
    VECTOR(vertex_recoding)[vertex]=1;
  }
  /* create vertex recoding vector */
  for (remaining_vertices=0, i=0; i<no_of_nodes; i++) {
    if (VECTOR(vertex_recoding)[i]==0) {
      VECTOR(vertex_recoding)[i]=remaining_vertices+1;
      remaining_vertices++;
    } else {
      VECTOR(vertex_recoding)[i]=0;
    }
  }
  /* create edge recoding vector */
  for (remaining_edges=0, i=0; i<no_of_edges; i++) {
    long int from=VECTOR(graph->from)[i];
    long int to=VECTOR(graph->to)[i];
    if (VECTOR(vertex_recoding)[from] != 0 &&
	VECTOR(vertex_recoding)[to  ] != 0) {
      VECTOR(edge_recoding)[i]=remaining_edges+1;
      remaining_edges++;
    } 
  }

  /* start creating the graph */
  newgraph.n=remaining_vertices;
  newgraph.directed=graph->directed;  

  /* allocate vectors */
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.from, remaining_edges);
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.to, remaining_edges);
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.oi, remaining_edges);
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.ii, remaining_edges);
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.os, remaining_vertices+1);
  IGRAPH_VECTOR_INIT_FINALLY(&newgraph.is, remaining_vertices+1);
  
  /* Add the edges */
  for (i=0, j=0; j<remaining_edges; i++) {
    if (VECTOR(edge_recoding)[i]>0) {
      long int from=VECTOR(graph->from)[i];
      long int to=VECTOR(graph->to  )[i];
      VECTOR(newgraph.from)[j]=VECTOR(vertex_recoding)[from]-1;
      VECTOR(newgraph.to  )[j]=VECTOR(vertex_recoding)[to]-1;
      j++;
    }
  }
  /* update oi & ii */
  IGRAPH_CHECK(igraph_vector_order(&newgraph.from, &newgraph.oi, 
				   remaining_vertices));
  IGRAPH_CHECK(igraph_vector_order(&newgraph.to, &newgraph.ii, 
				   remaining_vertices));  

  IGRAPH_CHECK(igraph_i_create_start(&newgraph.os, &newgraph.from, 
				     &newgraph.oi, remaining_vertices));
  IGRAPH_CHECK(igraph_i_create_start(&newgraph.is, &newgraph.to,
				     &newgraph.ii, remaining_vertices));
  
  /* attributes */
  IGRAPH_I_ATTRIBUTE_COPY(&newgraph, graph);
  IGRAPH_FINALLY_CLEAN(6);
  IGRAPH_FINALLY(igraph_destroy, &newgraph);
  IGRAPH_I_ATTRIBUTE_DELETE_VERTICES(&newgraph, &edge_recoding, 
				     &vertex_recoding);
  
  igraph_vit_destroy(&vit);
  igraph_vector_destroy(&edge_recoding);
  igraph_vector_destroy(&vertex_recoding);
  igraph_destroy(graph);
  *graph=newgraph;

  IGRAPH_FINALLY_CLEAN(4);
  return 0;
}