Warren Buffet says recession is going to be bad, invests in chewing gum!

"NEW YORK (Reuters) - Warren Buffett, the world's richest person, said on Monday the U.S. economy is in a recession that will be more severe than most people expect.

Buffett made his comments on CNBC television after his Berkshire Hathaway Inc (BRKa.N) (BRKb.N) agreed to invest $6.5 billion in the takeover of chewing gum maker Wm Wrigley Jr Co (WWY.N) by Mars Inc in a $23 billion transaction.

"This is not a field of specialty for me, but my general feeling is that the recession will be longer and deeper than most people think," Buffett said. "This will not be short and shallow.

"I think consumers are feeling gas and food prices," he added, "and not feeling they've got a lot of money for other things.""

http://news.yahoo.com/s/nm/20080428/bs_nm/buffett_recession_dc_1

Chewing gum prices are positively correlated with the prices of food - the more expensive the food is, the more chewing gum is used to replace it :-)...

The next development - George Soros shorting the weight-loss industry - to be announced later!

Computer Science lectures

This blog has a bunch of pointers to what looks like a nice collection of CS lectures.

http://freescienceonline.blogspot.com/2008/03/more-computing-video-lecture-courses.html

Election criteria

One of my favorite sayings: "Smart people tend to hire people smarter than themselves, dumb people tend to hire people dumber than themselves".

How does this work in the US? Let's see...

• We've got Bush for 2(!) terms.


• Elitist is a dirty word in politics.


• Incidentally, we've got Bush because Kerry was an elitist.


• The other reason we've got Bush is because a whole bunch of morons thought that Bush is more like them which is probably true.


• Hillary Clinton won Pennsylvania primary because voters without college education thought she is more like them. Which is both right and wrong - she's most certainly not like them, but she will definitely pander to them.


• There are no scientists in US Senate.


• But the most powerful Republican at the apex of Republican power in Congress was an exterminator. Yes: as in killing roaches.


• As a result, when confronted with statistics that clearly show that abstinence-only sex education in school does not work, a Republican congressman declares statistics "elitist".


• One third of US teenagers drops out of high school. Another third is not prepared for the next stage of life — either productive work or some form of post-secondary education".


• A quarter of US teenagers cannot identify Adolf Hitler. Large percentage of them think that US was allied with Germany against Russia during WWII. Both trends a boon to neocons because calling someone a "new Hitler" gets quick results from the dumbed down populace.


• Incidentally, there's probably not another place on earth where "neoCON" would be a legitimate name for a political current, not a stigma.


• Especially since they have actually perpetrated a classic Nigerian Scam leading to the Iraq War.

I think THE prerequisite for an effective democracy is an educated demos. This we definitely do not have.

US Aid to Israel

A few interesting numbers here: http://www.wrmea.com/html/us_aid_to_israel.htm

From Bill Maher

"If you think that Democrats are going to take your Bible away, you're an idiot.
If you think that they're going to take your gun away, you're an armed idiot.
And if you think that they're going to take your gun and give it to Mexicans to kill your god, you're Bill O'Reilly"

STL strings

So, while reviewing someone's code, I ran into a place where it was possible to eliminate an unnecessary string assignment. Which got me thinking - is mentioning this even worth it? How expensive is a string assignment, anyway?

So I wrote this simple program:

#define _SECURE_SCL 0

#include <string>

using namespace std;

void f(string a, string *b) {
  *b = a;
}

int main(int argc, char **argv) {
  string y = argv[0];
  string x;
  f(y, &x);
  return 0;
}

and stepped through it in the disassembler. The function f itself was inlined, so I only counted... prepare... the instructions in this call:

call        dword ptr [__imp_std::basic_string<char,std::char_traits<char>,std::allocator<char> >::operator= (402048h)] 

I skipped the call to new (argv[0] turns out to be 69 bytes, just above the default string buffer, so it does allocate), and memcpy itself (memcpy_s to be exact). The net result was - 270+ instructions. If there is no allocation, it is ~150!

This is correct - anywhere between 150 and 270 instructions of STL goo per string assignment IN ADDITION TO ACTUALLY DOING THE WORK!

I tried to do the same on Linux, and after 70-something ddd hung disassembling one of the functions...

Server 2008, first impressions

Short version

MUCH better than Vista.

Medium version

If you're a Microsoft employee and the goons from ITG(*) are trying to rip your favorite OS from your cold, dead hands and make you run Vista, don't - use Server 2008 instead.

There's no LUA/UAC/this idiotic thing that asks you whether you really intend to do what you just asked your computer to do. ipconfig works in a normal shell, instead of demanding one with elevated priorities. It takes literally half as many clicks/keystrokes to do everything.

There are manifestly fewer bugs. It did not fall apart within the first couple of weeks of use. And it seemed much snappier than Vista (although I did not do any real performance benchmarks).

Here is a blog post that has instructions on how to configure Server 2008 into a workstation: http://blogs.msdn.com/vijaysk/archive/2008/02/11/using-windows-server-2008-as-a-super-desktop-os.aspx..

And this one claims that it's 20% faster than Vista, but does not give details on how this number was arrived at: http://vista.blorge.com/2008/03/11/windows-server-2008-is-20-faster-than-vista/.
----
ITG - Microsoft's Informational Technology Group. Changed names multiple times over the years, but its essense stayed the same - these are the people who prevent developers from doing their jobs in the name of security. But this probably merits a separate post...

Long version

I have a lot of storage in my house - not counting client computers, there are approximately 9TB of redundant, RAIDed usable (after RAID) space. Most of my data is on this storage - software, music, pretty much everything.

Every time I buy anything, I immediately copy it to the server, and put the originals in a big box in the basement, where they will stay until the time when BSA people show up at my door and demand the proof of license :-).

My own data is replicated to multiple servers - approximately 300GB of home videos and digital photos, plus other, smaller stuff that accumulated over the years.

Storing stuff requires storage. Storing terabytes of things requires redundant storage. Over the years, I sampled a few RAID-5 solutions, both at home and at work.

There are two major problems with RAID controllers.
(1) While they protect you from a disk failure, they do not protect you from the failure of the controller itself.

Since all of these controllers use proprietary information to describe the RAID array that they store on the disks themselves, they are not interchangeable - you can't take a bunch of disks that you used in RAID mode on LSI and move them to Adaptec (while preserving that data that's on them, that is).

In fact, there's no guarantee that you can move disks between the controllers from the same manufacturer. Or between controllers with different versions of firmware. Or...

So 5 years from now when the RAID card fails, one can very easily be stuck with trying to find an exact replacement of a controller that had been out of production for the last 4.5 years. eBay, anyone?

(2) Software that accompanies these cards is often crappy.

The UI is almost always some atrocious Java program obviously written by a contractor in 2 days right before the product shipped, rife with misspellings and terrible English usage.

I have had multiple problems where midrange cards corrupted data when used on machines with more cores than the manufacturers originally expected.

So unless you test the disk failure scenario right upfront, BEFORE you get any data on it, you may well discover that (a) either the recovery mode does not work, or (2) because of unobvious UI you did something that wiped your disks instead of recovering them.

And good luck finding drivers when you upgrade to the new OS. And since you can't easily move the disks to a new controller... see (1).

Luckily, Microsoft Server family has software RAID subsystem (I am sure Linux has something similar, but coming from Microsoft, I am more familiar with Windows software).

To use it, you have to make your disks dynamic, then you can combine multiple disks into RAID-0, 1, or 5. Volumes of different types can share the same set of physical disks, so for example you can have part of disks 1 and 2, and whole of disk 3 to contribute to RAID-5 volume, and remainder of disks 1 and 2 to form RAID-0 temp storage.

The advantage of soft RAID is that it's hardware independent. You can take all these disks, shove them into any other computer (running Server), and it will still be a RAID volume with all your data intact. The same disk packs can be played on both Server 2003 and Server 2008.

It's an insanely cool idea. Unfortunately, in Server 2003, it was coupled with atrocious implementation.

Soft RAID-5 on Server 2003 is slow. Glacially slow. On my servers that feature dual Xeon 5130s (4 cores per server), with a very decent server motherboard (Supermicro X7DVL-E), and fairly decent mid-range SATA controllers, it was barely doing 20MBps writes, and sometimes would drop to 10MBps for extended periods of time. That on disks that are individually capable of 300MBps transfer rate.

RAID-5 works by partitioning disks into chunks, and then combining chunks from N disks to get N-1 chunk worth of data and 1 chunk of parity. Which means that to write a single sector to a volume, the RAID would have to read corresponding chunks from N disks, compute the parity, and write 2 sectors - one data and 1 parity.

A really terrible implementation would not cache the results of this read, so if the next sector needs to be written, it would repeat all the operations anew, instead of reusing the results of the previous reads.

The only way I can explain the RAID-5 write speed on Server 2003 is that it was this very terrible implementation, although I don't know for sure - an alternative explanation is that maybe they had sleep cycles in there :-).

So when Server 2008 came out, I could not wait to install it and check out its soft RAID implementation. I installed it first on my media server, and then on my data server.

Overall, I was quite impressed. Of course my expectations were very low to begin with because of Vista, but this thing was closer to Server 2003 than it was to Vista. I hit a few bugs right upfront - it hard hung once within a couple of days of installation, and then lost a set of disks (but recovered after reboot).

I am not quite ready to blame it on server itself though, because I added an unknown RAID controller to the machine, and it is more than likely that buggy drivers are to blame. The second server which did not have that controller did not (yet) exhibit this behavior.

Since then, it was relatively quiet and everything functioned the way it supposed to.

The drivers for SATA controllers from Server 2003 worked on Server 2008. The chipset drivers for the motherboard did not. I found that the chipset support for Server 2008 is still quite scarce.

What is unquestionably a bug in Server 2008 is that on RAID volumes the performance counters for logical disks are completely broken - everything except the disk free space and idle time is 0 when it is reading or writing full speed.

But most importantly, its soft RAID implementation is way faster than Server 2003. I copied a few terabytes of data so far, and on writes it does sustained throughput of ~80MBps - 4 times faster than the peak performance Server 2003 could muster. The reads (comparing files between two servers) almost saturate 1GBps network.

So far this things gets my stamp of approval :-).

I am yet to see if it is has long-term stability to last between Windows Update reboots (just in case, I preserved the original installations of Server 2003). I will report on this in a couple of months if everything goes well, earlier if it does not.

Henry Blodget agrees with me...

In today's post on Huffington Post he's basically saying the same thing I wrote about here: http://1-800-magic.blogspot.com/2007/12/risk-vs-reward.html: the CEO comp structure in the US encourages reckless, irresponsible attitude to the long-term shareholder value.

And he's very authoritative on the subject of irresponsible behavior :-)!

Who is flying this plane?

A while ago when I was in Technology Management MBA @ UW, they got former Qwest CEO Joseph Nacchio to come to the orientation meeting and talk about his tenure with the company.

As a result of this talk, I've gotten two very lasting impressions. First was a sense that under absolutely no conditions save outright starvation I would want to work for Qwest, and second was a state of bewilderment to which extent a CEO of the company can be disconnected from the technology that this company is using.

During his talk Nacchio confused GPRS, GSM, and other cell network terms. He professed the ignorance of how the networks work himself, and with a distinct sense of smugness (as in - you don't need to know any of this crap to be a leader - this is for lowly engineers to understand).

What he did talk about at length, and where he had become really animated and very lucid was mergers and acquisitions. Apparently, at the time Qwest was making a lot of money by buying small Eastern European telecoms and then selling them at a profit. From what I gathered in the talk, that line of business was what really animated the management. Not the pesky technological and operational hurdles of providing telecom services.

Over the last many years, as Ballmer's influence increased, I watched Microsoft leadership become less and less engineers, and more and more sales and marketing people.

Now, I don't think that engineering skills are required to head an engineering company. However, I do find that culturally marketing and engineering organizations are VERY different, almost stereotypically, Scott Adams, cartoonish different.

And nothing affects company's culture more than the CEO.

So when I ran into this "pearl" today, I felt embarrassed, but, sadly, not surprised:



There is often a big disconnect between market's perception of the product, and the sales perception it. It is not bad - to be an effective salesperson, you have to be excited about the product as it is - or you won't sell very many of it. Sadly, I fear that Microsoft internal perception - that of the leadership team, anyway - of Vista is closer to the movie above, and not to the painful reality.

Google blog post about the incidents

Here's the post on the Google Lat Long blog about the traffic incidents:
http://google-latlong.blogspot.com/2008/04/new-ways-to-beat-traffic-with-google.html.

Google Road Traffic Incidents ahoy!

I have not posted for a long, long while because the confluence of several events has been keeping me busy almost around the clock for the last 5 weeks.

First and foremost, the project I've been working for the last 3 months - the road traffic incidents - has shipped today. It's my first experience building the entire pipeline at Google - downloading the data from 3rd party provider, parsing it, storing the results in a bigtable, writing a spatial index for it, and, finally, publishing it through the maps frontend server. Here's the result:


It's actually kinda fun watching the data in various cities. Seattle is a really, really quiet place. The incidents are mostly construction (there are a few accident icons, but in reality they are all about slow traffic):


Though Seattle is a small place, actual real accidents (as in - collisions) are quite uncommon even in much bigger cities. Here's New York. Here you see a lot of road closures (most of them are periodic, and in the future - eventually we'll figure out how to be more intelligent about displaying them; right now the data we're getting is lacking a lot of details that would allow us to be - but we'll work it out with our provider). There are a few constructions, and slowdowns, but only one of the icons is displaying a collision:


And most other US cities are somewhere between New York and Seattle. Except for one. If you're living in LA, my hat is off to you. The traffic there is real zoo. Yes, there's even one place where almost everyday the sign shows "Animals on the road". Here's a quiet night at LA (believe it or not, during the day time, it's a lot worse). Most non-consturction icons are real collisions. A couple are hit and runs.


The insane amount of construction icons you see here is a nightly event - the authorities deply miriads of crews to do various road maintenance tasks on the roads in the evening, and it lasts through the night. And they report a construction incident for every one of them. During the day, the roadwork in LA clears, leaving mostly just collisions. But a lot of them!

Adam Smith on wars...

"In great empires the people who live in the capital, and in the provinces remote from the scene of action, feel, many of them, scarce any inconveniency from the war; but enjoy, at their ease, the amusement of reading in the newspapers the exploits of their own fleets and armies. To them this amusement compensates the small difference between the taxes which they pay on account of the war, and those which they had been accustomed to pay in time of peace. They are commonly dissatisfied with the return of peace, which puts an end to their amusement, and to a thousand visionary hopes of conquest and national glory from a longer continuance of the war."

-- An Inquiry into the Nature And Causes of the Wealth of Nations

Dr. Strangelove... in reverse!

President Bush rarely comments about the Democratic presidential contest, but he said that he had to speak up about Clinton's red phone ads because he found them "so confusing."

"If I answered the red phone every time it rang, I would never get any sleep," Bush said. "Sometimes it starts ringing at 9 p.m., and I am already tucked in by then."

Bush said that "there's nothing so important that it can't wait until tomorrow, or whenever I remember to check my voicemail."


http://news.yahoo.com/s/uc/20080308/cm_uc_crabox/op_475477;_ylt=AmW7uqcqRkkMlQ.hJ9SIda9xFb8C

This is a loop that never ends - it just goes on and on my friends!

    for (double d = 0.0; d != 1.0; d += 0.1)
      System.out.println(d);

:-)

10,000 BC: a movie review

Short version: terrible.

Medium version: think 300. This movie is very similar to it in terms of both connection to reality and the overall flow.

Long version (WARNING: a plot spoiler): there's this really advanced civilization in the middle of desert 12000 years ago. They know metalworks (and it looks like iron, too), shipbuilding, navigation, astronomy, writing, cloth manufacture. Basically, Atlantis, except in Africa. They build huge gold-tipped pyramids, elevated roads, and babylonian-looking temples and are rulled by a priest class reporting to chief priest AKA the "Almighty". To build all this they need, obviously, slaves, and so they send the roaming bands to snatch them from neighboring tribes.

The neighbors are more authentically looking 10000BC-ers (at least they wear fur, and appear to not use iron weapons). The only problem - their men are either shaved, or have neatly trimmed beards and goatees.

Anyways, one of these roving bands steals this good-looking girl, and it proves to be their undoing: her boyfriend comes to rescue her, incites the global slave revolt, easily overcomes the few armed (iron swords!) guards that are there. In the thick of the fight the gilfriend gets shot, but then is magically revived by a magician thousands of miles away.

Oh, yes, and at the end of all this our hero's tribe gets a bag of magic beans and starts the agricultural revolution. Viva Monsanto!

Java collections perf redux

I finally got some time to run my microbenchmark for various Java collection objects.

Here's the source:

package com.solyanik.perftesting;

import java.util.Arrays;
import java.util.Vector;
import java.util.ArrayList;

/**
 * Implements a micro-benchmark using
 * high-resolution RDTSC timer.
 * 
 * @author sergey@solyanik.com (Sergey Solyanik)
 */
public class MicroBenchmark {

  /*
   * RDTSC - returns number of clocks since
   * the core was last reset.
   * 
   * NOTE: no synchronizing instruction
   * is executed before RDTSC, so this
   * is not suitable for really short
   * runs as the second RDTSC can execute
   * before the instruction we're trying
   * to measure.
   * 
   * NOTE2: Different cores have
   * different clock counters, and they
   * are not synchronized.
   */
  private native static long rdtsc();
  
  /**
   * A function that consists of one
   * instruction - ret.
   */
  private native static void nop();
  
  /**
   * Fixes the thread to the first CPU
   * and core.
   */
  private native static void fixthread();

  /**
   * Loads the library containing our
   * native code.
   */
  static {
    System.loadLibrary("jnipc");
  }
  
  /**
   * Prints the statistics for an array
   * of clock intervals.
   * 
   * @param clocks The array of clocks
   * @param iterations Number of iterations
   * @param bucket_number How many buckets
   *     to use for the histogram
   * @param scale Scale for the histogram
   */
  private static void Analyze(
      String header, long [] clocks,
      int iterations, int bucket_number,
      int scale) {
    Arrays.sort(clocks);
    double average = 0.0;

    long [] buckets = new long[bucket_number];

    double log2 = Math.log(2.0);

    for (int i = 0 ; i < clocks.length; ++i) {
      average += clocks[i];
      int ndx = (int)(Math.log(
          (double)clocks[i]
          / (double)scale) / log2);
      if (ndx < 0)
        buckets[0]++;
      else if (ndx < buckets.length - 1)
        buckets[ndx + 1]++;
      else
        buckets[buckets.length - 1]++;
    }

    average /= clocks.length;

    System.out.println(header + ":");
    System.out.println("Median for " +
        iterations + " iterations: " +
        clocks[clocks.length / 2]);
    System.out.println("Average for " +
        iterations + " iterations: " +
        average);
    System.out.println("Median per iteration: "
        + clocks[clocks.length / 2]
                 / iterations);
    System.out.println("Average per iteration: "
        + average / iterations);

    System.out.println("Histogram:");    
    int prev = 0;
    int curr = scale;
    for (int i = 0; i < buckets.length - 1; ++i) {
      System.out.println(prev + " - " + curr
          + ": " + buckets[i]);
      prev = curr;
      curr *= 2;
    }

    System.out.println(" > " + prev + ": " +
        buckets[buckets.length - 1]);    
  }
  
  /**
   * Number of loop iterations per sample.
   */
  private static final int ITERATIONS = 10000;
  
  /**
   * Number of samples.
   */
  private static final int SAMPLES = 1000;
  
  /**
   * Number of buckets for the histogram.
   */
  private static final int BUCKETS = 20;
  
  /**
   * Bucket scale for the histogram.
   */
  private static final int BUCKET_SCALE = 1000;
  
  /**
   * Main function.
   * 
   * @param args
   */
  public static void main(String[] args) {
    fixthread();
    System.out.println("Let the fun begin!");
    
    long [] clocks = new long[SAMPLES];

    Boolean b = new Boolean(true);

    for (int i = 0; i < clocks.length; ++i) {
      Boolean [] a = new Boolean[ITERATIONS];
      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        a[j] = b;
      long stop = rdtsc();
      clocks[i] = stop - start;
    }
    
    Analyze("Array",
        clocks, ITERATIONS, BUCKETS,
        BUCKET_SCALE);
        
    for (int i = 0; i < clocks.length; ++i) {
      Vector<Object> v = new Vector<Object>();

      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        v.add(b);
      long stop = rdtsc();
      clocks[i] = stop - start;
    }

    Analyze("Expanding Vector", clocks,
        ITERATIONS, BUCKETS, BUCKET_SCALE);

    for (int i = 0; i < clocks.length; ++i) {
      Vector<Object> v = new Vector<Object>();
      v.setSize(ITERATIONS);
      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        v.set(j, b);
      long stop = rdtsc();
      clocks[i] = stop - start;
    }

    Analyze("Pre-allocated Vector",
        clocks, ITERATIONS, BUCKETS,
        BUCKET_SCALE);

    for (int i = 0; i < clocks.length; ++i) {
      ArrayList<Object> al = new ArrayList<Object>();
      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        al.add(b);
      long stop = rdtsc();
      clocks[i] = stop - start;
    }

    Analyze("Expanding ArrayList", clocks,
        ITERATIONS, BUCKETS, BUCKET_SCALE);

    for (int i = 0; i < clocks.length; ++i) {
      ArrayList<Object> al = new ArrayList<Object>();
      al.ensureCapacity(ITERATIONS);
      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        al.add(b);
      long stop = rdtsc();
      clocks[i] = stop - start;
    }

    Analyze("Pre-allocated ArrayList", clocks,
        ITERATIONS, BUCKETS, BUCKET_SCALE);
  }
}

And here's what I've got:

                     Median Average
Array                    14      15
Expanding vector        124     133
Pre-allocated vector    115     119
Expanding ArrayList      57      66
Pre-allocated ArrayList  43      46

To remind, C++ numbers were like this:

                     Median Average
Array                     2       2
Pre-allocated vector,
  assigment               2       2
Pre-allocated vector,
  push_back               8       8
Expanding vector,
  push_back             145     148 (*)

(*) Obviously, this number very much depends on number of expansions, because STL vector grows by 50%, a very big vectors will have smaller cost. See my previous posts for details on C++ testing.

Synchronization in a web application

UPDATED 03/08/2008 - source code completely rewritten.

While I was writing the code for the Project Guttenberg proxy, one of the problems I needed to solve was the concurrency in URL processing.

When a user requests a URL from the proxy, the proxy sends the request to the real Project Gutenberg web site, and writes the result into a file.

Then it processes the file - if it's an HTML document, it replaces the references to the original web site with links that point to the proxy. If it was a request for a LRF file, the text file that was fetched from the Project Gutenberg would be translated.

Then there's cache - the web pages that have already been fetched, but have not expired get fetched from that cache.

In other words, there's a lot of activity, that is by and large independent - unless the two requests come for the same page. Then of course we do not want to process them at the same time, because they are going to be modifying the same state, corrupting each other in the process.

The solution is not that hard, but I thought it would be instructive to put it here.

package com.solyanik.sync;

import java.util.HashSet;

/** The general object lock class.
 * Ensures that any two objects
 * that are equal to each other as
 * in a.equals(b) cannot be locked
 * simultaneously.
 * 
 * This is a complement to the
 * standard monitor functionality
 * that protects objects that are
 * equal in == sense (i. e. the
 * same object).
 *
 * IMPORTANT NOTE: The same thread
 * trying to lock the same object
 * twice = INSTANT DEADLOCK!
 * 
 * Usage:
 *    // in constructor
 *    ObjectLock lock = new ObjectLock();
 *    ... 
 *    void handleRequest(final Object o) {
 *      try {
 *        lock.lock(o);
 *        ... // do processing
 *      } finally {
 *        lock.unlock(o);
 *      }
 *
 * @author sergey@solyanik.com (Sergey Solyanik)
 */
class ObjectLock {
  /**
   * Stores all currently locked objects
   */
  private HashSet<Object> locked =
      new HashSet<Object>();
  
  public void lock(Object o)
      throws InterruptedException {
    synchronized(locked) {
      for ( ; ; ) {
        if (locked.contains(o)) {
          // We are already locked.
          locked.wait();
        } else {
          locked.add(o);
          break;
        }
      }
    }
  }
  
  public void unlock(Object o)
      throws IllegalMonitorStateException {
    synchronized(locked) {
      if (!locked.remove(o))
        throw new
            IllegalMonitorStateException(
            "The object is not locked.");
      locked.notifyAll();
    }
  }
}

The code is tiny and by and large self explanatory. We keep a hash set, and whenever we are about to lock the object, we check if it is in the hash set already. If it is not, we put it there and continue on our merry way.

If it is, we're joining the queue of threads waiting on the hash set (it can actually be absolutely any object). When the object is unlocked, it will be removed from the hash, and the threads waiting will be notified.

To preempt obvious question - this assumes that the level of contention is not super high, so there are not a lot of threads pending at any given time. In the case of a high level of contention (10s of threads), we would want to wake up only the threads that are waiting for that specific object that has just been unlocked, not all of the threads.

How is this different from synchronized(o)? The ObjectLock will prevent collisions of code operating on objects that are equal in the sense that o1.equals(o2). The synchronized clause prevents operations on the objects that are the same (o1 == o2).

In the case of a string, by the way, an inefficient way of doing the same that the lock class above is trying to accomplish is this:

  void handle_request(String url) {
    String lock = url.toLowerCase().intern();
    synchronized(lock) {
       ... do processing ...
    }
  }

Interned string have the property of s1 == s2 iff s1.equals(s2). The problem of course is that every url that has ever been processed will stay in your process's string heap forever, so it is not very efficient.

And the unit test:

package com.solyanik.sync;

import java.util.HashSet;

public class UnitTest {
  static class Executor extends Thread {
    private static HashSet<String> inuse
        = new HashSet<String> ();

    private String url;
    private int id;
    private ObjectLock lock;
    
    public Executor(String url, int id,
        ObjectLock lock) {
      this.url = url;
      this.id = id;
      this.lock = lock;
    }

    public void run() {
      try {
        System.out.println("[" + id +
            "] Started to process " + url);
        System.out.println("[" + id +
            "] Locking " + url);
        lock.lock(url);
        System.out.println("[" + id +
            "] Lock acquired " + url);
        // Yes, I know this is unnecessary.
        // But just to be absolutely sure,
        // this is a test after
        // all...
        String internedUrl = url.intern();

        if (inuse.contains(internedUrl)) {
          System.out.println("[" + id +
              "] ERROR: COLLISION " + url);
          throw new AssertionError(
              "Bad lock implementation!");
        }

        inuse.add(internedUrl);
        Thread.sleep(5000);
        inuse.remove(internedUrl);

        System.out.println("[" + id +
            "] Releasing " + url);
        lock.unlock(url);
        System.out.println("[" + id +
            "] Lock released " + url);
        System.out.println("[" + id +
            "] Thread finished ");
      } catch (InterruptedException ex) {
        System.out.println("[" + id +
            "] Interrupted: " + ex);
      }
    }
  }
  
  /**
   * @param args Not used.
   */
  public static void main(String[] args) {
    String [] urls = {
       new String("http://www.google.com"),
       new String("http://www.google.com"),
       new String("http://www.google.com"),
       new String("http://www.google.com"),
       new String("http://www.google.com"),
       new String("http://www.live.com"),
       new String("http://www.live.com"),
       new String("http://www.live.com"),
       new String("http://www.yahoo.com"),
       new String("http://www.yahoo.com")
    };
    ObjectLock lock = new ObjectLock();
    for (int i = 0; i < 128; ++i) {
      (new Executor(urls[i % urls.length],
          i, lock)).start();
    }
  }
}

JNI Made Simple

I want to measure the performance of Java vector to compare it with STL. Java of course lacks a really high resolution timer, so I needed to use JNI - the Java native interface.

After reading Sun's spec here: http://java.sun.com/j2se/1.5.0/docs/guide/jni/spec/jniTOC.html, I decided that there needs to be a simplified version. So here it goes. A simple, 3-step recipe.

(1) Create a DLL. Here's an example of a DLL that contains one function that does nothing:

nop.cxx:

__declspec(naked)
unsigned __int64 __cdecl nop(void) {
  __asm {
    ret
  }
}

nop.def:
LIBRARY "nop"
EXPORTS
  Java_Test_nop = nop

Notice the decorated name. Java_ is required for all, Test is the name of the class, "nop" is the name of the function. Straightforward.

(2) Compile this using Visual Studio, and copy the resulting nop.dll to the root directory of your Java project.

(3) Call it from Java:

class Test {
  native static void nop();
  
  static {
    System.loadLibrary("nop");
  }

  public static void main(String[] args) {
    nop();
  }
}

That's all there is to it. Let's now look at how fast the function call is actually made. To do this, I am going to use the same framework I used for testing STL's vector performance: http://1-800-magic.blogspot.com/2008/02/stl-vector-performance.html, except translated to Java.

We'll need 3 functions in the native library though - one that does nothing (this is the function we're going to test), the other will be returning the results of Pentium's RDTSC instruction that measures the clock count, and the third one will be fixing the thread's affinity to the first core, because if this is not done and thread jumps the core while executing, we won't be able to compare the results of RDTSC, since clock counters on multiple cores is not synchronized.

Here's how the native code looks like:

jnipc.cxx:
#include <windows.h>

__declspec(naked)
unsigned __int64 __cdecl rdtsc(void) {
  __asm {
    rdtsc
    ret
  }
}

__declspec(naked)
unsigned __int64 __cdecl nop(void) {
  __asm {
    ret
  }
}

void __cdecl fixthread(void) {
  SetThreadAffinityMask(GetCurrentThread(), 1);
}


jnipc.def:
LIBRARY "jnipc"
EXPORTS
  Java_Test_rdtsc = rdtsc
  Java_Test_nop = nop
  Java_Test_fixthread = fixthread

The Java code is more involved, but less complicated :-):

import java.util.Arrays;

class Test {
  native static long rdtsc();
  native static void nop();
  native static void fixthread();
  
  static {
    System.loadLibrary("jnipc");
  }
  
  static final int ITERATIONS = 100000;
  static final int SAMPLES = 1000;
  static final int BUCKETS = 20;
  static final int BUCKET_SCALE = 10000;
  
  public static void main(String[] args) {
    fixthread();
    long [] clocks = new long[SAMPLES];
    for (int i = 0; i < SAMPLES; ++i) {
      long start = rdtsc();
      for (int j = 0; j < ITERATIONS; ++j)
        nop();
      long stop = rdtsc();
      clocks[i] = stop - start;
    }

    Arrays.sort(clocks);
    double average = 0.0;

    long [] buckets = new long[BUCKETS];

    double log2 = Math.log(2.0);

    for (int i = 0 ; i < SAMPLES; ++i) {
      average += clocks[i];
      int ndx = (int)(Math.log((double)clocks[i]
          / (double)BUCKET_SCALE) / log2);
      if (ndx < 0)
        buckets[0]++;
      else if (ndx < BUCKETS-1)
        buckets[ndx + 1]++;
      else
        buckets[BUCKETS-1]++;
    }

    average /= SAMPLES;

    System.out.println("Median for " +
        ITERATIONS + " iterations: " +
        clocks[SAMPLES / 2]);
    System.out.println("Average for " +
        ITERATIONS + " iterations: " +
        average);
    System.out.println("Median per iteration: "
        + clocks[SAMPLES / 2] / ITERATIONS);
    System.out.println("Average per iteration: "
        + average / ITERATIONS);

    System.out.println("Histogram:");    
    int prev = 0;
    int curr = BUCKET_SCALE;
    for (int i = 0; i < BUCKETS - 1; ++i) {
      System.out.println(prev + " - " + curr
          + ": " + buckets[i]);
      prev = curr;
      curr *= 2;
    }

    System.out.println(" > " + prev + ": " +
        buckets[BUCKETS - 1]);

  }
}

On my T2600 Core 2 Duo laptop this ran in approximately 38 clocks per iteration.

Median for 100000 iterations: 3815903
Average for 100000 iterations: 3854832.709
Median per iteration: 38
Average per iteration: 38.54832709
Histogram:
0 - 10000: 0
10000 - 20000: 0
20000 - 40000: 0
40000 - 80000: 0
80000 - 160000: 0
160000 - 320000: 0
320000 - 640000: 0
640000 - 1280000: 0
1280000 - 2560000: 0
2560000 - 5120000: 998
5120000 - 10240000: 2
10240000 - 20480000: 0
20480000 - 40960000: 0
40960000 - 81920000: 0
81920000 - 163840000: 0
163840000 - 327680000: 0
327680000 - 655360000: 0
655360000 - 1310720000: 0
1310720000 - -1673527296: 0
 > -1673527296: 0


C:\bin>java -version
java version "1.6.0_03"
Java(TM) SE Runtime Environment (build 1.6.0_03-b05)
Java HotSpot(TM) Client VM (build 1.6.0_03-b05, mixed mode, sharing)

To compare it with C, I've added the following case to my perf program in http://1-800-magic.blogspot.com/2008/02/stl-vector-performance.html:

case 4:
  printf("Function call\n");
  HMODULE hm = LoadLibraryW(L"c:\\bin\\jnipc.dll");
  void (*func)() = (void (*)())
    GetProcAddress(hm, "Java_Test_nop");
  for (int i = 0; i < NUM_SAMPLES; ++i) {
    unsigned __int64 begin = rdtsc();
    for (int j = 0; j < NUM_ITERATIONS; ++j)
      func();
    unsigned __int64 end = rdtsc();
    clocks[cptr++] =
        (unsigned int)(end - begin);
}

...and got approximately 6 clocks per call, making JNI about 6 times slower on this simple test:

Function call
Median: 60112
Median, CPI: 6
Average: 60262
Average, CPI: 6
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 0
32000 - 64000: 998
64000 - 128000: 1
128000 - 256000: 1
256000 - 512000: 0
512000 - 1024000: 0
1024000 - 2048000: 0
2048000 - 4096000: 0
4096000 - 8192000: 0
8192000 - 16384000: 0
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0

Testing disk backup

Now that I've written a simple disk backup problem, I need to ensure that it is correct. After all, we're dealing with the user data here - terrible things will happen if it will read (or write) a wrong cluster once in a while!

So, obviously, before I can use it, I need to persuade myself that it won't eat my data.

What makes a great test? The best test of all runs the software on all possible data sets in all possible environments, and verifies that it ran correctly. In our case, verification means that the backup can be restored, and the restored files match the files that were backed up perfectly.

Verification is actually is quite straightforward - we'll take a backup of one partition, restore it to another, and then compare the files.

Much harder is to generate all combinations of file systems and all environments (or partition sizes, in this case).

When there's a test that needs a lot of various combination, and I can't foresee the "interesting" data sets, I like the random number generator. I like it for two reasons. First, given enough time, it will come up with things I could not foresee. Second, most random number generators aren't really random - they generate a fixed sequence of numbers, given the starting number (a seed). For testing, this is actually an asset - this means that the test that is driven by a pseudorandom number generator is REPRODUCIBLE. Why is this good (actually, necessary)? Because when the test does not work, you want to repeat the conditions under which the software broke.

What about running in all possible environments? For us, the environment is the partition size. Here we will have to yield to reality - there is no way we can test a lot of partition configurations. We'll just have to go with the sequence I consider "interesting" - partition sizes of powers of two, powers of two +1, powers of two -1, and a random size in between.

Ok, then, the test. It was written in python and run from the same directory where the source for the program is.

"""Test for discio.

Usage:
  python.exe test.py disk_no letter1 letter2
Where:
  disk_no The number of the disk to use.
      The disk must contain no information.
      Everything on the disk will be wiped out!
  letter1 The letter for the first volume to be
      created on the test disk. Must not be
      an existing volume, or else it will be
      wiped out.
  letter2 The letter for the second volume to be
      created on the test disk. Must not be
      an existing volume, or else it will be
      wiped out.

  WARNING: If run incorrectly, this test can
      be enormously destructive! It cleans disks,
      formats volumes etc. If you call it with an
      incorrect disk number, that disk will be gone!
"""

import os
import random
import sys
import tempfile
import time

class Error(Exception):
  """Base error for the test."""
  pass


class GenericFailure(Error):
  """Test has failed. The string describes the failure."""
  pass


__DISKPART = 'c:\\windows\\system32\\diskpart.exe'
__FORMAT = 'c:\\windows\\system32\\format.com'
__DISKIO = 'Release\\discio.exe'
__BACKUPFILE = 'test.bkf'
__SVI = 'System Volume Information'

A function that compares directories and files. If given two files, it just keeps reading them in 2K chunks, and comparing the results. If given two directories, it lists them, sorts them, and then compares the sorted lists in the merge-sort fashion: if the names are the same, match the files (or directories). If not, say that the smaller (in lexicographical order) does not exists, and the advance the pointer to the smaller one.

def Compare(f1, f2):
  """Compares two filesystem objects.

  Prints out the differences.
  Returns:
    True if the files compare OK.
  """
  ret = True

  if (__SVI in f1 and
      __SVI in f2):
    return ret

  if os.path.isfile(f1) and os.path.isfile(f2):
    fd1 = open(f1, 'rb')
    fd2 = open(f2, 'rb')
    while True:
      s1 = fd1.read(2048)
      s2 = fd2.read(2048)
      if s1 != s2:
         ret = False
         break
      if not s1:
         break
    fd1.close()
    fd2.close()
    if not ret:
      print '%s and %s are different' % (f1, f2)
    return ret

  if os.path.isfile(f2):
    print '%s is a dir, %s is a file' % (f1, f2)
    return False

  if os.path.isfile(f1):
    print '%s is a dir, %s is a file' % (f2, f1)
    return False

  #print "%s <-> %s" % (f1, f2)

  files1 = os.listdir(f1)
  files2 = os.listdir(f2)

  files1.sort()
  files2.sort()

  i = 0
  j = 0

  while i < len(files1) and j < len(files2):
    if files1[i] == files2[j]:
      if not Compare(os.path.join(f1, files1[i]),
                     os.path.join(f2, files2[j])):
        ret = False
      i = i + 1
      j = j + 1
      continue

    ret = False

    if files1[i] < files2[j]:
      print 'Missing ' + os.path.join(f2, files1[i])
      i = i + 1
    else:
      print 'Missing ' + os.path.join(f1, files2[j])
      j = j + 1

  while i < len(files1):
    print 'Missing ' + os.path.join(f2, files1[i])
    i = i + 1
    ret = False

  while j < len(files2):
    print 'Missing ' + os.path.join(f1, files2[j])
    j = j + 1
    ret = False

  return ret    

This one is the biggest bottleneck in the whole system: imagine generating 4G of data and writing it to files in a SCRIPT! Yet that's exactly what it does, computer time be damned!

def GenerateFiles(base, size):
  """Generates a bunch of random files and directories.

  Args:
    base The base directory where the files should be.
        Must end with '\\'.
    size The size in kilobytes of files that it should
        not exceed.
  """
  if not os.path.exists(base):
    os.mkdir(base)

  if size > 100:
    GenerateFiles(os.path.join(base, 'dir1'), size / 6)
    GenerateFiles(os.path.join(base, 'dir2'), size / 6)
    GenerateFiles(os.path.join(base, 'dir3'), size / 6)
    size = size / 2

  file_num = 0
  size = size * 1024

  while size > 0:
    if size < 2000:
      file_size = size
    else:
      file_size = random.randint(0, size)
    fd = open(os.path.join(base, ('file%d.txt' % file_num)),
              'w')

    str = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
    str = str + str.lower() + '\n'

    try:
      for i in xrange(0, file_size/len(str)):
        fd.write(str)

      fd.flush()
    finally:
      fd.close()

    file_num = file_num + 1
    size = size - file_size

The reason I remove some files is to generate holes - unused clusters between the used ones. Otherwise I would be backing up a completely filled disk - almost as uninsteresting as testing backup of an empty disk!

def RemoveSomeFiles(base):
  """Removes every third file from a tree."""
  counter = 0
  for (basepath, dirs, files) in os.walk(base):
    if __SVI in basepath:
      continue
    for file in files:
      if counter % 3 == 0:
        full_file_name = os.path.join(basepath, file)
        os.remove(full_file_name)
      counter = counter + 1


def GenerateFileSystem(base, size):
  """Generates a bunch of random files and directories.

  Args:
    base The base directory where the files should be.
        Must end with '\\'.
    size The size in kilobytes of files that it should
        not exceed.
  """
  try:
    GenerateFiles(base, int(size * 9 / 10))
  except IOError:
    pass

  RemoveSomeFiles(base)

Finally, the test itself. We use diskpart to wipe and repartition the disk on every iteration - which means that we need the whole disk to ourselves to run. If you haven't figured this out on your own by now :-).

def RunTest(size, fs, disk, source, dest):
  """The test.

  Creates a two partitions, on the disk,
  formats one, fills it with a bunch of random
  files, backups it, and restores the results
  on the second partition.

  Then verifies that the files are identical
  on both.

  Args:
    size The size of an individual partition
    fs The file system to use, FAT, FAT32,
        or NTFS
    disk The number of the disk to use in testing
        Everything on this disk will be gone!
    source The letter to use for volume on which
        the test data will be generated
    dest The letter to use for volume to which
        the backup will be restored.

  Throws:
    GenericFailure if the test fails

  """
  print 'Running the test.'
  print '  Physical disk: %d ' % disk
  print '  Parition size: %dM' % size
  print '  File system:   %s' % fs

  (fi, script) = tempfile.mkstemp()
  os.close(fi)
  fd = open(script, 'w')
  fd.write('SELE DIS %d\n' % disk)
  fd.write('CLEAN\n')
  fd.write('CREAT PART PRI SIZE=%d\n' % size)
  fd.write('ASSIGN LETTER=%s\n' % source)
  fd.write('CREAT PART PRI SIZE=%d\n' % size)
  fd.write('ASSIGN LETTER=%s\n' % dest)
  fd.write('EXIT\n')
  fd.close()

  print 'Preparing the partitions'

  ret = os.spawnv(os.P_WAIT, __DISKPART,
                  [__DISKPART,
                   '/s',
                   script])

  os.remove(script)

  if ret != 0:
    raise GenericFailure(
        'Could not partition the disk!')

  ret = os.spawnv(os.P_WAIT, __FORMAT,
                  [__FORMAT,
                   source + ':',
                   '/fs:' + fs,
                   '/v:test',
                   '/y'])

  if ret != 0:
    raise GenericFailure(
        'Could not format the test partition')

  print 'Generating test files'

  GenerateFileSystem(source + ':\\', size * 1024)

  time.sleep(20)

  print 'Running backup'

  ret = os.spawnv(os.P_WAIT, __DISKIO,
                  [__DISKIO,
                   source + ':\\',
                   'readto',
                   __BACKUPFILE]);

  if ret != 0:
    raise GenericFailure('Backup failed!')

  print 'Running restore'

  ret = os.spawnv(os.P_WAIT, __DISKIO,
                  [__DISKIO,
                   dest + ':\\',
                   'writefrom',
                   __BACKUPFILE]);

  if ret != 0:
    raise GenericFailure('Restore failed!')

  if not Compare(source + ':\\', dest + ':\\'):
    raise GenericFailure(
        'Restore does not match backup!')

  print 'Test succeeded!'
  print

And finally, the main. First and foremost it is trying to keep you safe: prevent you from wiping out some useful data.

def main(args):
  """Main entry point.

  Args:
    args: Command line arguments. See usage.

  """

  random.seed(0)

  if len(args) < 3:
    print __doc__
    return
    
  try:
    disk = int(args[0])
  except ValueError:
    print 'Disk should be a number!'
    return

  if disk <= 0:
    print 'Disk cannot be a system disk (0).'
    return

  source = args[1].upper()
  dest = args[2].upper()
  if len(source) != 1 or len(dest) != 1:
    print 'Only use 1 letter to specify drives'
    return

  if (dest < 'K' or dest > 'O'
      or source < 'K' or source > 'O'):
    print 'Use letters K through O for test drives.'
    print 'This is to minimize potential of a'
    print 'collision with the real drives.'
    return

  if dest == source:
    print 'Source and destination should be different!'
    return

  if os.path.exists(source + ':\\'):
    print 'Cannot continue: %s:\\ exists!' % source
    return

  if os.path.exists(dest + ':\\'):
    print 'Cannot continue: %s:\\ exists!' % dest
    return

  print 'Running the test on physical disk %d' % disk
  print 'and using drive letters %s and %s' % (source,
                                               dest)

And here we run through a bunch of partition sizes (in Megabytes) I thought were interesting, and run backup test on 3 file systems: FAT16, FAT32, and NTFS.

  try:
    for size in [50, 100, 128, 200, 256, 300, 315,
                 1000, 1024, 2047, 2048, 2049, 4095,
                 4096, 4097]:
      RunTest(size, "FAT", disk, source, dest)
      RunTest(size, "FAT32", disk, source, dest)
      RunTest(size, "NTFS", disk, source, dest)
  except GenericFailure, f:
    print "Test failed: " + f.message

  
if __name__ == '__main__':
  main(sys.argv[1:])

And - voila! This ran overnight and succeeded!

Taking disk image snapshots - programmatically

A while ago I published a quick program to rip software CD images by opening CD device as a file for block access.

Today I will expand it to read (and restore) the images of a regular hard drive.

Technically speaking, just running the code from http://1-800-magic.blogspot.com/2008/02/backing-up-software-cds.html directly would work (sort of, we'll talk about an important gotcha later). But it will result in a very, very big file. Hard drives are huge these days, but mostly empty. So a lot of the backup image would be the empty space.

Luckily, Windows has a handy facility to figure out what blocks on the disks are actually used. Device IO control FSCTL_GET_VOLUME_BITMAP returns a bitmap where 1s correspond to used blocks, and 0s to the clusters that are empty.

Unluckily, this io control only works for NTFS. For FAT it returns the bitmap, but it does not start with the beginning of the disk - rather, with the end of the FAT tables as it seems.

Funnily enough, FSCTL_GET_VOLUME_BITMAP also reports the starting cluster where the bitmap is supposed to start. For FAT, it should have been non-zero - we could then backup everything before it, and then whatever is in the bitmap after, but no luck - it just happily returns zero, and the shifted bitmap.

There's surely a way to figure out the sizes of FATs, but since this filesystem is less and less common, and the drives formatted with FAT are small, so in the interests of simplicity I just backup the whole disk if the filesystem is not NTFS.

Another bit of interesting data is that your partition is actually bigger than your file system by one cluster. This cluster contains a copy of the the boot sector for the partition. The intent, I think, is to use it in disk recovery. We back it up for NTFS only, not for other file systems.

Alright, so here's the code. The usual things first. See this post for implementation of AtoI: http://1-800-magic.blogspot.com/2008/02/down-with-atoi.html.

#define UNICODE 1
#define _UNICODE 1

#include <windows.h>
#include <stdio.h>
#include <ctype.h>
#include <strsafe.h>
#include <assert.h>

#define ARRLEN(c) (sizeof(c)/sizeof(c[0]))

template<class IntType> bool AtoI(const WCHAR *sz, IntType *out) {
...}

GetDeviceName translates the disk number (0-N), a drive letter (A:), or a volume mount point (C:\, or c:\cdrom) to a name that Windows would understand as a block device.

bool GetDeviceName(WCHAR *szDevName,
                   size_t cchDevName,
                   const WCHAR *szInput) {
  unsigned int disk_no = 0;
  if (AtoI<unsigned int>(szInput, &disk_no)) {
    // It's a physical disk
    StringCchPrintfW(szDevName,
                     cchDevName,
                     L"\\\\.\\PHYSICALDRIVE%d", 
                     disk_no);
    return true;
  }

  WCHAR drive_letter = toupper(szInput[0]);
  if (drive_letter >= 'A' && 
      drive_letter <= 'Z' &&
      szInput[1] == ':' &&
      szInput[2] == '\0') {
    // It's a drive letter
    StringCchPrintfW(szDevName,
                     cchDevName,
                     L"\\\\.\\%c:",
                     drive_letter);
    return true;
  }

  WCHAR sz[_MAX_PATH];
  const WCHAR *p = wcsrchr(szInput, '\\');
  if (!p || p[1] != '\0') {
    // Mount point needs to end in backslash
    StringCchCopyW(sz, ARRLEN(sz), szInput);
    StringCchCatW(sz, ARRLEN(sz), L"\\");
    szInput = sz;
  }

  if (!GetVolumeNameForVolumeMountPointW(
      szInput, szDevName, cchDevName))
    return false;

  // This returns a name ending with '\'
  // but CreateFile does not take it, so
  // we need to get rid of it.
  WCHAR *q = wcsrchr(szDevName, '\\');
  if (q && q[1] == '\0')
    *q = '\0';

  return true;
}

Then opening is straightforward. Notice the sharing mode - this is part of the gotcha I will be talking about later. Also notice that I removed the FSCTL_ALLOW_EXTENDED_DASD_IO from here. This IO control allows reading and writing past the end of the volume as understood by the file system. Since we backup the overflow sector only for NTFS, I moved it to NTFS-specific place.

HANDLE OpenBlockDevice(const WCHAR *szDevName,
                       bool fWrite) {
  HANDLE h = CreateFile(szDevName,
    fWrite ? GENERIC_WRITE : GENERIC_READ,
    FILE_SHARE_WRITE, NULL, OPEN_EXISTING,
    FILE_FLAG_NO_BUFFERING |
        (fWrite ? FILE_FLAG_WRITE_THROUGH : 0),
    NULL);
  return h;
}

Now the fun part - let's get the volume bitmap. The idiosyncrasy of this function is that there's no way - at least not a way I know of - to get the size of the bitmap that we need to pre-allocate, first. If you give it the NULL input, as is customary with majority of Windows API, and expect the output to say how many bytes it needs, you will be disappointed - the function just errors out, and returns 0 in the bytes returned. Therefore, I do this ugly dance:

VOLUME_BITMAP_BUFFER *GetVolumeBitmap(
    HANDLE hDev, unsigned int *puiBitmapSize) {
  STARTING_LCN_INPUT_BUFFER sStartLcn;
  sStartLcn.StartingLcn.QuadPart = 0;

  DWORD dwBitmapSize = 0;
  DWORD dwAllocatedSize = 64 * 1024;
  VOLUME_BITMAP_BUFFER *pVolumeBitmap = NULL;
  for ( ; ; ) {
    pVolumeBitmap = (VOLUME_BITMAP_BUFFER *)
        LocalAlloc(LMEM_FIXED, dwAllocatedSize);

    BOOL ret = DeviceIoControl(hDev,
                      FSCTL_GET_VOLUME_BITMAP,
                      &sStartLcn,
                      sizeof(sStartLcn),
                      pVolumeBitmap,
                      dwAllocatedSize,
                      &dwBitmapSize,
                      NULL);
    if (ret) {
      *puiBitmapSize = dwBitmapSize;
      return pVolumeBitmap;
    }

    if (GetLastError() != ERROR_MORE_DATA)
      return NULL;

    dwAllocatedSize *= 2;
  }
}

I want to support both the case where we backup the entire volume as well as the case where we would only backup the used clusters (let's call it sparse backup). I do it by prepending the following header to the sparse backup file - so when I read it back, if the magic string is there, I know that I need to read the volume bitmap, and then write the data accordingly. If it's not there, the file is basically the image and I just blast the whole thing on the disk:. When I restore, I also need the sizes of the cluster and the sector, so I know how to interpret the bitmap. If I were somewhat less lazy, I would also care to check that on restore the sector size of the disk matches that of the backup.

#define SIG "Compressed Disk Image"

struct BackupHeader {
  char signature[sizeof(SIG)];
  unsigned int uiSectorSize;
  unsigned int uiClusterSize;
  unsigned int uiBitmapSize;
};

I put together a generic function that can transfer data from disk to file, and from file to disk according to the volume bitmap if it is given. The small gotcha here is that the size of the volume bitmap is in BITS. Storage API rarely follows standard Windows patters for some reason!

Notice that in the case where pVolumeBitmap is NULL, the function gets reduced to a very simple loop that just reads and writes the entire block device.

// This function assumes the ownership of
// pVolumeBitmap, i. e. it frees it.
// hFrom is the handle from which the data is read.
// hTo is the handle to which the data is written.
// hSeek is always the handle of the device - this
// is the pointer that is being moved to skip empty
// extents.
// For backup: hFrom = hSeek - disk handle
//             hTo - file handle
// For restore: hTo = hSeek - disk handle
//             hFrom - file handle
// If device is not an NTFS volume, pBackupData,
// pVolumeBitmap are both NULL.
DWORD Transfer(HANDLE hFrom,
               HANDLE hTo,
               HANDLE hSeek,
               BackupHeader *pBackupData,
               VOLUME_BITMAP_BUFFER *pVolumeBitmap,
               void *buffer,
               const unsigned int cbBufferSize,
               unsigned __int64 &ui64Bytes) {
  DWORD dwErr = ERROR_SUCCESS;

  // This is used to step through the bitmap
  // of used clusters.
  unsigned char bmpMask = 1;
  unsigned char *bmpIndex = NULL;
  int maxClustersPerRead = 0;
  unsigned __int64 maxClusters = 0;
  if (pVolumeBitmap) {
    bmpIndex = pVolumeBitmap->Buffer;
    maxClusters = pVolumeBitmap->BitmapSize.QuadPart;
    maxClustersPerRead = cbBufferSize /
        pBackupData->uiClusterSize;
  }

  unsigned __int64 cluster = 0;
  bool fTerm = false;

  while (!fTerm) {
    DWORD dwToMove = cbBufferSize;
    if (pVolumeBitmap) {
      // Skip empty. Note that at the end
      // of this loop cluster can be more
      // than maxCluster, because sometimes
      // we step in 8.
      while (cluster < maxClusters) {
        if (!*bmpIndex) {
          cluster += 8;

          ++bmpIndex;
          assert(bmpMask == 1);

          continue;
        }

        if ((*bmpIndex & bmpMask) == 0) {
          ++cluster;

          bmpMask <<= 1;
          if (! bmpMask) {
            bmpMask = 1;
            ++bmpIndex;
          }

          continue;
        }
        break;
      }

      dwToMove = 0;
      int numClusters = 0;
      while (cluster  + numClusters < maxClusters
          && numClusters < maxClustersPerRead
          && (*bmpIndex & bmpMask) != 0) {
        dwToMove += pBackupData->uiClusterSize;
        ++numClusters;

        bmpMask <<= 1;
        if (!bmpMask) {
          bmpMask = 1;
          ++bmpIndex;
        }
      }

      assert(dwToMove <= cbBufferSize);
      if (dwToMove == 0) { // The End?
        // The last sector on the volume
        // contains duplicate of the boot
        // sector, but is not part of the
        // volume map.
        dwToMove = pBackupData->uiSectorSize;

        // Since we might have skipped more
        // clusters than there were because
        // sometimes we count by 8-s, we reset
        // the current cluster count as well:
        cluster = maxClusters;

        // Convenient place to free the bitmap.
        LocalFree(pVolumeBitmap);
        pVolumeBitmap = NULL;
        fTerm = true;
      }

      unsigned __int64 offset = cluster *
        (unsigned __int64)pBackupData->uiClusterSize;

      if (!SetFilePointerEx(hSeek,
                            *(LARGE_INTEGER *)&offset,
                            NULL,
                            FILE_BEGIN)) {
        dwErr = GetLastError();
        wprintf(L"Seek error %d (0x%08x)\n",
            dwErr, dwErr);
      }

      cluster += numClusters;
    }

    DWORD dwRead = 0;
    if (!ReadFile(hFrom, buffer,
        dwToMove, &dwRead, NULL)) {
        dwErr = GetLastError();
        wprintf(L"\nRead error %d (0x%08x)\n",
            dwErr, dwErr);
        break;
    }

    if (dwRead == 0)
      break;

    DWORD dwWrit = 0;
    if (!WriteFile(hTo, buffer,
        dwRead, &dwWrit, NULL) ||
        dwWrit != dwRead) {
      dwErr = GetLastError();
      wprintf(L"\nWrite error %d (0x%08x)\n",
          dwErr, dwErr);
      break;
    }

    ui64Bytes += dwWrit;
    wprintf(L"%I64u\r", ui64Bytes);
  }

  return dwErr;
}

This makes read much shorter than what it would have been:

DWORD ReadToFile(HANDLE hDev,
                 HANDLE hFile,
                 void *buffer,
                 const unsigned int cbBufferSize,
                 unsigned __int64 &ui64BytesRead) {
  DWORD dwErr = ERROR_SUCCESS;

  unsigned int uiBitmapSize = 0;
  VOLUME_BITMAP_BUFFER *pVolumeBitmap = NULL;

  NTFS_VOLUME_DATA_BUFFER sVolumeData;
  DWORD dwRead = 0;
  // For NTFS volumes we can store just the
  // clusters that are used. We get the
  // bitmap where 1 means that the cluster
  // is used here.
  if (DeviceIoControl(hDev,
                      FSCTL_GET_NTFS_VOLUME_DATA,
                      NULL,
                      0,
                      &sVolumeData,
                      sizeof(sVolumeData),
                      &dwRead,
                      NULL)) {
    // For NTFS disable volume boundary checks by FS
    DWORD dwBytes = 0;
    DeviceIoControl(hDev,
                    FSCTL_ALLOW_EXTENDED_DASD_IO,
                    NULL,
                    0,
                    NULL,
                    0,
                    &dwBytes,
                    NULL);

    pVolumeBitmap = GetVolumeBitmap(hDev, &uiBitmapSize);
  }

  BackupHeader h;
  DWORD dwWrit = 0;
  if (pVolumeBitmap) {
    memcpy(h.signature, SIG, sizeof(SIG));
    h.uiBitmapSize = uiBitmapSize;
    h.uiSectorSize = sVolumeData.BytesPerSector;
    h.uiClusterSize = sVolumeData.BytesPerCluster;
    if (!WriteFile(hFile,
                   &h,
                   sizeof(h),
                   &dwWrit,
                   NULL)
        || dwWrit != sizeof(h)) {
      dwErr = GetLastError();
      wprintf(L"Write error %d (0x%08x)\n",
          dwErr, dwErr);
      return dwErr;
    }

    if (!WriteFile(hFile,
                   pVolumeBitmap,
                   uiBitmapSize,
                   &dwWrit,
                   NULL)
        || dwWrit != uiBitmapSize) {
      dwErr = GetLastError();
      wprintf(L"Write error %d (0x%08x)\n",
          dwErr, dwErr);
      return dwErr;
    }
  }

  return Transfer(hDev,
                  hFile,
                  hDev,
                  &h,
                  pVolumeBitmap,
                  buffer,
                  cbBufferSize,
                  ui64BytesRead);
}

And write, even shorter:

DWORD WriteFromFile(HANDLE hDev,
                 HANDLE hFile,
                 void *buffer,
                 const unsigned int cbBufferSize,
                 unsigned __int64 &ui64BytesWritten) {
  DWORD dwRead = 0;
  DeviceIoControl(hDev,
                  FSCTL_ALLOW_EXTENDED_DASD_IO,
                  NULL,
                  0,
                  NULL,
                  0,
                  &dwRead,
                  NULL);

  unsigned int uiBitmapSize = 0;
  VOLUME_BITMAP_BUFFER *pVolumeBitmap = NULL;

  BackupHeader h;

  if (ReadFile(hFile, &h, sizeof(h), &dwRead, NULL)
      && dwRead == sizeof(h)) {
    if (memcmp(h.signature, SIG, sizeof(SIG)) == 0) {
      pVolumeBitmap = (VOLUME_BITMAP_BUFFER *)
          LocalAlloc(LMEM_FIXED, h.uiBitmapSize);
      if (!ReadFile(hFile, pVolumeBitmap, h.uiBitmapSize,
          &dwRead, NULL) || dwRead != h.uiBitmapSize) {
        wprintf(L"Backup corrupted - could not read the "
                L"volume bitmap!\n");
        return ERROR_DATA_NOT_ACCEPTED;
      }
    } else
      SetFilePointer(hFile, 0, NULL, FILE_BEGIN);
  }

  DWORD dwErr = Transfer(hFile,
                         hDev,
                         hDev,
                         &h,
                         pVolumeBitmap,
                         buffer,
                         cbBufferSize,
                         ui64BytesWritten);

  DeviceIoControl(hDev,
                  IOCTL_DISK_UPDATE_PROPERTIES,
                  NULL,
                  0,
                  NULL,
                  0,
                  &dwRead,
                  NULL);

  return dwErr;
}

The main just parses the command line then:

int wmain(int argc, WCHAR **argv) {
  if (argc < 4) {
    wprintf(L"Usage: %s {disk|volume} cmd "
            L"args\n", argv[0]);
    wprintf(L"Where disk is physical number "
            L"of the drive, e. g. 0\n");
    wprintf(L"Volume is a drive letter or a "
            L"full path to \n");
    wprintf(L"the mount point, e. g. a: or "
            L"c:\\mount\\vol\n");
    wprintf(L"Commands:\n");
    wprintf(L"readto filename - read the "
            L"contents of the device into "
            L"a file\n");
    wprintf(L"writefrom filename - write the "
            L"contents of the file onto "
            L"the disk or volume\n");
    return 0;
  }

  UINT uiErrModeSav = SetErrorMode(
      SEM_FAILCRITICALERRORS);

  WCHAR szDevName[_MAX_PATH];
  if (!GetDeviceName(szDevName,
                     ARRLEN(szDevName),
                     argv[1])) {
    wprintf(L"Could not recognize %s\n", argv[1]);
    SetErrorMode(uiErrModeSav);
    return 1;
  }

  const int kBufferSize = 64 * 1024;
  void *buffer = VirtualAlloc(
      NULL,
      kBufferSize,
      MEM_COMMIT | MEM_RESERVE,
      PAGE_READWRITE);
  if (!buffer) {
    DWORD dwErr = GetLastError();
    wprintf (L"Could not reserve memory, "
             L"error %d (0x%08x)\n", dwErr, dwErr);
    SetErrorMode(uiErrModeSav);
    return 3;
  }

  DWORD dwTickStart = GetTickCount();
  unsigned __int64 ui64Bytes = 0;
  DWORD dwErr = ERROR_SUCCESS;
  if (_wcsicmp(argv[2], L"readto") == 0) {
    HANDLE hDev = OpenBlockDevice(szDevName, false);
    if (hDev != INVALID_HANDLE_VALUE) {
      HANDLE hFile = CreateFile(
          argv[3], GENERIC_WRITE, 0, NULL,
          CREATE_ALWAYS, 0, NULL);
      if (hFile != INVALID_HANDLE_VALUE) {
        wprintf(L"%s --> %s:\n", szDevName, argv[3]);
        dwErr = ReadToFile(hDev,
                           hFile,
                           buffer,
                           kBufferSize,
                           ui64Bytes);
        CloseHandle(hFile);
        if (dwErr != ERROR_SUCCESS)
          DeleteFile(argv[3]);
      } else {
        dwErr = GetLastError();
        wprintf (L"Could not open %s, error %d "
                 L"(0x%08x)\n", argv[3],
                 dwErr, dwErr);
      }
      CloseHandle(hDev);
    } else {
      dwErr = GetLastError();
      wprintf(L"Could not open %s, "
              L"error %d (0x%08x)\n",
              argv[1], dwErr, dwErr);
    }
  } else if (_wcsicmp(argv[2], L"writefrom") == 0) {
    HANDLE hDev = OpenBlockDevice(szDevName, true);
    if (hDev != INVALID_HANDLE_VALUE) {
      HANDLE hFile = CreateFile(
          argv[3], GENERIC_READ, 0, NULL,
          OPEN_EXISTING, 0, NULL);
      if (hFile != INVALID_HANDLE_VALUE) {
        wprintf(L"%s --> %s:\n", argv[3], szDevName);
        dwErr = WriteFromFile(hDev,
                              hFile,
                              buffer,
                              kBufferSize,
                              ui64Bytes);
        CloseHandle(hFile);
      } else {
        dwErr = GetLastError();
        wprintf (L"Could not open %s, error %d "
                 L"(0x%08x)\n", argv[3],
                 dwErr, dwErr);
      }
      CloseHandle(hDev);
    } else {
      dwErr = GetLastError();
      wprintf(L"Could not open %s, "
              L"error %d (0x%08x)\n",
              argv[1], dwErr, dwErr);
    }
  } else {
    wprintf (L"Command %s not recognized.\n",
             argv[2]);
    dwErr = ERROR_NOT_SUPPORTED;
  }

  if (dwErr == ERROR_SUCCESS &&
      ui64Bytes != 0) {
    wprintf(L"Transferred %I64u bytes in %d "
            L"seconds.\n", ui64Bytes,
            (GetTickCount() - dwTickStart) / 1000);
  }

  VirtualFree(buffer, 0, MEM_RELEASE);
  SetErrorMode(uiErrModeSav);

  return dwErr == ERROR_SUCCESS ? 0 : 4;
}

Now it WOULD be time to publish the binary and the source code in its intact form, except for one thing. This code works poorly if it is run against your system volume. The problem is that the system volume does not usually exist in the consistent state - it is being constantly written to, with part of the information already on disk, and another part in buffers.

Moreover, while the backup is in progress, another program can add data, thus creating an inconsistency between the disk and the NTFS descriptors (say, we'd read the disk, but not yet the MFT. Then the program would add a file to the area of the disk we have already written, and modify the MFT - we would've gotten an image where there's a file full of garbage).

The solution is called VSS, and I will cover it in another post. Hopefully, at that point we will have the code that would work well with the system volume. But even before that, I owe an article on testing it.

As it stands now, it can be used on secondary volumes (when you're sure nothing is writing on it), unprotected CDs and DVDs, or on a system volume if you're booting an alternative OS.

STL vector perf redux

Alright, my exercise in STL vector perf analysis was to a large extent an exercise in n00bness. Nevertheless, I learned a lot of stuff, thanks to the readers, and here's the summary.

If you are looking for the gory details, read this: http://1-800-magic.blogspot.com/2008/02/stl-vector-performance.html and the commentary.

The task: stuff an STL vector with 10000 numbers, measure what it takes per iteration.

The performance of course will vary with the computer. Mine is T2600 @ 2.16GHz.

C reference point - pre-allocated array.
2 clocks per iteration:

int *work = (int *)malloc(10000 * sizeof(int));
for (int j = 0; j < 10000; ++j)
  work[j] = j;

C reference point - dynamically allocated array using STL allocation policy.
70 clocks per iteration:

int size = 20;
int *work = (int *)malloc(size * sizeof(int));
for (int j = 0; j < 10000; ++j) {
  if (j >= size) {
    size = size + size / 2;
    work = (int *)realloc(work,
        size * sizeof(int));
  }
  work[j] = j;
}

STL, dynamically allocated vector, with or without security checks.
140 clocks per iteration.

vector work;
for (int j = 0; j < 10000; ++j)
  work.push_back(j);

STL, pre-allocated, with security checks, push_back.
9 clocks per iteration.

vector work;
work.reserve(10000);
for (int j = 0; j < 10000; ++j)
  work.push_back(j);

STL, pre-allocated, no security checks, push_back.
8 clocks per iteration.

#define _SECURE_SCL 0
vector work;
work.reserve(10000);
for (int j = 0; j < 10000; ++j)
  work.push_back(j);

STL, pre-allocated, with security checks, [].
6 clocks per iteration.

vector work(10000);
for (int j = 0; j < 10000; ++j)
  work[j] = j;

STL, pre-allocated, no security checks, [].
2 clocks per iteration.

#define _SECURE_SCL 0
vector work(10000);
for (int j = 0; j < 10000; ++j)
  work[j] = j;

Resume: know what you're doing. Performance can be from 3x to 70x slower if you're not careful.

For those who didn't read the previous post, and want to argue that none of this is important, read this http://1-800-magic.blogspot.com/2007/12/memory-is-not-free-more-on-vista.html for a fun real life story when you wouldn't think that performance matters until the very end, when it suddenly does, but then it's too late :-).

Another observation: the fact that most STL implementations are not THAT much smaller is an accolade to the modern CPU design.

If you actually look at the generated code, the STL (especially push_back) compiles to much bigger assembly than you could guess from the results above - even without the security check it is ~30 instructions plus the call into reallocating logic.

; for (int j = 0; j < NUM_ITERATIONS; ++j)
mov         edx,dword ptr [esp+3Ch] 
mov         ecx,dword ptr [esp+38h] 
xor         esi,esi 
mov         dword ptr [esp+1Ch],esi 
mov         edi,edi 
; work.push_back(j);
L0:
cmp         ecx,edi 
jne         L1
xor         eax,eax 
jmp         L2
L1:
mov         eax,dword ptr [esp+40h] 
sub         eax,ecx 
sar         eax,2 
L2:
mov         ebx,edx 
sub         ebx,ecx 
sar         ebx,2 
cmp         ebx,eax 
jae         L3
mov         dword ptr [edx],esi 
add         edx,4 
mov         dword ptr [esp+3Ch],edx 
jmp         L4
L3:
lea         ecx,[esp+1Ch] 
push        ecx  
push        edx  
lea         eax,[esp+3Ch] 
call        std::vector<int,std::allocator<int> >::_Insert_n
mov         edx,dword ptr [esp+3Ch] 
mov         ecx,dword ptr [esp+38h] 
L4:
inc         esi  
cmp         esi,2710h 
mov         dword ptr [esp+1Ch],esi 
jl          $LN22+5Ah (4013E0h)

This is only 4 times slower than the C loop, which is so much less:

xor         eax,eax 
lea         esp,[esp]
L0: 
mov         dword ptr [ecx+eax*4],eax 
inc         eax  
cmp         eax,2710h 
jl          L0

This is because all these heaps of STL code execute in parallel, speculatively, whereas the straightforward C loop is not parallelizeable at all.

STL vector performance

I was going through an interview package today, where the interviewee returned a whole STL vector as a result of the function. I wondered if this is a reasonable thing to do, and whether it would cause a lot of data copying.

So I cranked up Visual Studio 2008 jotted down a few lines of code, and looked at the disassembly.

Turns out, the compiler does a really good job with the aliasing, so the code that it generates (in optimized version) is very efficient. What drew my attention was nearly 30 instructions that it generated for "vector.push_back(i)". I expected that the optimized code would be a lot less, and started digging.

As it turned out, it was a lot worse. The code generated inline was the tip of an iceberg - another function was called, and not just in the case where array reallocation was required. I stepped through it in assembly several times to make sure.

It looked pretty bad, so I decided to measure how it would compare to non-STL programming idioms. I compared pushing 10000 numbers into pre-allocated STL vector, an auto-expanding STL vector, an expanding C array, and a pre-allocated C array.

The code was compiled in Release (optimized) mode on Visual Studio 2008.

Here is the test. Getting the boring stuff out of the way first:

#include <windows.h>

#include <stdio.h>
#include <math.h>

#include <vector>

using namespace std;

I used Pentium clock counter instruction to minimize overhead. This function returns 64-bit clock count since the CPU was last reset. Thank you, Wikipedia!

__declspec(naked)
unsigned __int64 __cdecl rdtsc(void) {
  __asm {
    rdtsc
    ret
  }
}

However, if the thread is relocated to a different core, the clock count returned is for a different CPU, which makes time measurement completely random. So I restricted the code to a single core, and for a good measure, pre-committed the heap to ensure that paging does not spoil the results:

#define NUM_ITERATIONS 10000
void PrepEnvironment(void) {
  // Move everything to one CPU, we're going
  // to be using its clock counter
  SetThreadAffinityMask(GetCurrentThread(), 1);

  // Commit all the memory.
  void *p = malloc(sizeof(int) * 4 * NUM_ITERATIONS);
  memset(p, 0, sizeof(int) * 4 * NUM_ITERATIONS);
  free(p);
}

For stats, I calculated the average, the median, and the histogram on the log scale:

int __cdecl comparator(const void *p1,
                       const void *p2) {
  const unsigned int *pi1 =
      (const unsigned int*)p1;
  const unsigned int *pi2 =
      (const unsigned int*)p2;
  return (int)(*pi1 - *pi2);
}

#define BUCKETS 20
void DumpStats(unsigned int *clocks,
               unsigned int cptr,
               unsigned int iters) {
  qsort(clocks, cptr,
      sizeof(unsigned int), comparator);

  int buckets[BUCKETS];
  memset(buckets, 0, sizeof(buckets));

  double average = 0;
  double log2 = log(2.0);

  for (unsigned int i = 0 ; i < cptr; ++i) {
    average += (double)clocks[i];
    int ndx = (int)(log((double)clocks[i]
        / 1000.0) / log2);
    if (ndx < 0)
      buckets[0]++;
    else if (ndx < BUCKETS-1)
      buckets[ndx + 1]++;
    else
      buckets[BUCKETS-1]++;
  }

  average /= cptr;

  printf("Median: %d\n",
      (int)clocks[cptr / 2]);
  printf("Median, CPI: %d\n",
      clocks[cptr / 2] / iters);
  printf("Average: %d\n", (int)average);
  printf("Average, CPI: %d\n",
      (int)(average / iters));

  printf("Histogram:\n");
  int prev = 0;
  int curr = 1000;
  for (int i = 0; i < BUCKETS - 1; ++i) {
    printf("%d - %d: %d\n", prev, curr,
        buckets[i]);
    prev = curr;
    curr *= 2;
  }

  printf(" > %d: %d\n", prev,
      buckets[BUCKETS - 1]);
}

Finally, the 4 different pieces of code all filled an array with 10000 numbers, and repeated this operation 1000 times (so the standard error is about 3%). I have made C array reallocation strategy to match STL vector implementation.

#define NUM_SAMPLES 1000
int main(int argc, char **argv) {
  PrepEnvironment();

  static unsigned int clocks[NUM_SAMPLES];
  for (int k = 0; k < 4; ++k) {
    unsigned int cptr = 0;

    switch(k) {
    case 0:
      printf("STL re-alloc\n");
      for (int i = 0; i < NUM_SAMPLES; ++i) {
        vector work;
        unsigned __int64 begin = rdtsc();
        for (int j = 0; j < NUM_ITERATIONS; ++j)
          work.push_back(i);
        unsigned __int64 end = rdtsc();
        clocks[cptr++] =
            (unsigned int)(end - begin);
      }
      break;
    case 1:
      printf("STL pre-alloc\n");
      for (int i = 0; i < NUM_SAMPLES; ++i) {
        vector work;
        work.reserve(NUM_ITERATIONS);
        unsigned __int64 begin = rdtsc();
        for (int j = 0; j < NUM_ITERATIONS; ++j)
          work.push_back(i);
        unsigned __int64 end = rdtsc();
        clocks[cptr++] =
            (unsigned int)(end - begin);
      }
      break;
    case 2:
      printf("C re-alloc\n");
      for (int i = 0; i < NUM_SAMPLES; ++i) {
        int size = 20;
        int *work = (int *)malloc(size * sizeof(int));
        unsigned __int64 begin = rdtsc();
        for (int j = 0; j < NUM_ITERATIONS; ++j) {
          if (j >= size) {
            // Match STL vector allocation algorithm
            size = size + size / 2;
            work = (int *)realloc(work, size
                * sizeof(int));
          }
          work[j] = j;
        }
        unsigned __int64 end = rdtsc();
        clocks[cptr++] =
            (unsigned int)(end - begin);
        free(work);
      }
      break;
    case 3:
      printf("C pre-alloc\n");
      for (int i = 0; i < NUM_SAMPLES; ++i) {
        int *work = (int *)malloc(
            NUM_ITERATIONS * sizeof(int));
        unsigned __int64 begin = rdtsc();
        for (int j = 0; j < NUM_ITERATIONS; ++j)
          work[j] = j;
        unsigned __int64 end = rdtsc();
        clocks[cptr++] =
            (unsigned int)(end - begin);
        free(work);
      }
      break;
    }

    DumpStats(clocks, cptr, NUM_ITERATIONS);
    printf("\n\n");
  }

  return 0;
}

So the results are as follows (this is on my T60p ThinkPad T2600 Core 2 Duo running at 2.16 GHz plugged-in):

STL re-alloc
Median: 1511328
Median, CPI: 151
Average: 1540473
Average, CPI: 154
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 0
32000 - 64000: 0
64000 - 128000: 0
128000 - 256000: 0
256000 - 512000: 0
512000 - 1024000: 123
1024000 - 2048000: 835
2048000 - 4096000: 35
4096000 - 8192000: 6
8192000 - 16384000: 1
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0


STL pre-alloc
Median: 80171
Median, CPI: 8
Average: 82007
Average, CPI: 8
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 0
32000 - 64000: 0
64000 - 128000: 995
128000 - 256000: 2
256000 - 512000: 1
512000 - 1024000: 2
1024000 - 2048000: 0
2048000 - 4096000: 0
4096000 - 8192000: 0
8192000 - 16384000: 0
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0


C re-alloc
Median: 698529
Median, CPI: 69
Average: 728341
Average, CPI: 72
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 0
32000 - 64000: 0
64000 - 128000: 0
128000 - 256000: 0
256000 - 512000: 0
512000 - 1024000: 972
1024000 - 2048000: 27
2048000 - 4096000: 1
4096000 - 8192000: 0
8192000 - 16384000: 0
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0


C pre-alloc
Median: 20124
Median, CPI: 2
Average: 21572
Average, CPI: 2
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 980
32000 - 64000: 15
64000 - 128000: 2
128000 - 256000: 2
256000 - 512000: 0
512000 - 1024000: 1
1024000 - 2048000: 0
2048000 - 4096000: 0
4096000 - 8192000: 0
8192000 - 16384000: 0
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0

The results are quite damning: the most commonly accepted C idiom of pre-allocating the array is 2 clocks per integer inserted, compared to 8 for pre-allocated vector, and 150 for most commonly used dynamically expanding vector. Which is 4 to 75 times(!) faster.

NOTE: The original version of the article contained an error: I pre-allocated the array incorrectly. The problem is now fixed and results updated. The fix resulted in dramatic improvement of STL performance, although it still did not come close to matching raw C performance.

NOTE2: I have tested it with the disabled security (#define _SECURE_SCL 0, and verified that generated assembly did not contain calls to the boundary checker), but it did not change the performance.

Even the dynamically expanding C array is the same speed as the pre-allocated vector, and 2.5 times faster than its dynamically-expanding STL counterpart.

But who cares, right? Computers are fast, we'll optimize the tight loops, and nobody would notice things that are not on the critical path, anyway. Well, yes and no. If your software is slow across the board, you'll be hit with what I call the Vista syndrome - there isn't one loop to optimize if it's all sluggish. More on this here: http://1-800-magic.blogspot.com/2007/12/memory-is-not-free-more-on-vista.html.

For a company like Google, it also means higher CPU utilization, more machines, and ultimately, more energy used. The difference in power by the way is very real - Google clusters do have a certain maximum CPU utilization threshold, beyond which the fuses get blown.

Next on my plate is to measure Java performance, and my gut feel is that it's not going to be much slower than STL. Which brings up the question - why use C++ at all, if the only advantage it has over Java and C# - the performance - is eaten by the runtime?


-----
Update: There IS a way to get STL to match C in the pre-allocated case. Here's how it's done. Ouch!

...
#define _SECURE_SCL 0
#include 
...

  printf("STL pre-alloc\n");
  for (int i = 0; i < NUM_SAMPLES; ++i) {
    vector work(NUM_ITERATIONS);
    work.reserve(NUM_ITERATIONS);
    unsigned __int64 begin = rdtsc();
    for (int j = 0; j < NUM_ITERATIONS; ++j)
      work[j] = j;
    unsigned __int64 end = rdtsc();
    clocks[cptr++] =
        (unsigned int)(end - begin);
  }

...and here are the results:

STL pre-alloc
Median: 20124
Median, CPI: 2
Average: 21668
Average, CPI: 2
Histogram:
0 - 1000: 0
1000 - 2000: 0
2000 - 4000: 0
4000 - 8000: 0
8000 - 16000: 0
16000 - 32000: 979
32000 - 64000: 11
64000 - 128000: 5
128000 - 256000: 4
256000 - 512000: 1
512000 - 1024000: 0
1024000 - 2048000: 0
2048000 - 4096000: 0
4096000 - 8192000: 0
8192000 - 16384000: 0
16384000 - 32768000: 0
32768000 - 65536000: 0
65536000 - 131072000: 0
131072000 - 262144000: 0
 > 262144000: 0