Posts tagged with fortran

FORTRAN 语言是怎么来的

在高级语言是怎么来的子系列的第一篇中, 我们结合当时硬件的特点,分析了 FORTRAN 为什么一开始不支持递归。但是 FORTRAN 本身是怎么来的这个问题其实还是没有得到正面回答,本节我们就谈谈 FORTRAN 语言本身是怎么来的。

其实,FORTRAN 语言也是现实驱动的。 所以我们还是回到当时,看看当时程序员的需求和软硬件条件,看看 FORTRAN 是怎么来的。 了解历史的另一个好处是, 因为 FORTRAN 的发展历史正好和高级语言的发展历史高度重合,所以了解 FORTRAN 的背景,对于理解其他高级语言的产生都是大有帮助的。

1. 困难的浮点计算
我们先从硬件的角度说起。 大致从 1946 年第一台计算机诞生,到 1953 年,计算机一直都缺少两件非常重要的功能,一个叫浮点计算,一个叫数组下标寻址,这两个功能的缺失直接导致了高级语言的兴起。 我们依次单个分析。读者对浮点计算应该都不陌生,用通俗的话说就是如 0.98×12.6 这样的实数乘法,或者  0.98 + 12.6 这样的实数加法的运算。用行话说,就是用计算机进行大范围高精度数的算术运算。

学过二进制的同学都知道,二进制整数之间的乘法和加法等运算都是相对简单的,和正常的十进制运算是一样的,只是把加法和乘法这些基本操作用更加简单的逻辑或(OR) 和 逻辑与 (AND) 实现而已,在电子电路上也很好实现。 因此,就是世界上最早的电子计算机,ENIAC,也是支持整数的乘法加法等算术操作的。

可是浮点运算就不一样了。 因为一个额外的小数点的引入,在任何时候都要注意小数点的对齐。 如果用定点计数,则计数的范围受到限制,不能表示非常大或者非常小的数。所以,浮点数一般都是用科学记数法表示的,比如 IEEE 754 标准。(不熟悉 IEEE 754 的读者也可以想像一下如何设计一套高效的存储和操作浮点数的规范和标准,以及浮点算法), 科学记数法表示的浮点数的加减法每次都要对齐小数点,乘除法为了保持精度,在设计算法上也有很多技巧,所以说,相比较于整数的运算和逻辑运算,浮点运算是一件复杂的事情。落实到硬件上就是说,在硬件上设计一个浮点运算,需要复杂的电路和大量的电子元器件。在早期电子管计算机中,是很少能做到这么大的集成度的。因此,不支持浮点也是自然的设计取舍。在计算机上放一个浮点模块这个想法,需要等电子工业继续发展,使得电子管体积小一点,功耗低一点后,才能进入实践。

2. 关于浮点计算的一些八卦

关于浮点,这里顺带八卦一点浮点计算的事情。在计算机芯片设计中,浮点计算一直是一个让硬件工程师头疼的事情,即使到了386时代,386 处理器 (CPU)的浮点乘法也是用软件模拟的,如果想用硬件做浮点乘法,需要额外购买一块 80387 浮点协处理器 FPU,否则就在 386 上做软件的模拟。这样做的原因在一块硅片上刻蚀一个 CPU 和一个FPU 需要的集成度还是太高,当时的工艺根本没法实现。真的把 FPU 和 CPU 融在一起刻蚀到一块硅片上,已经是 1989 年的事情了。当时,Intel 把融合了 80386 和 80387 的芯片改了改,起了个名字叫 80486,推向了市场。带着浮点的处理器的普及,使得个人计算机能做的事情变多了。极度依赖于浮点计算的多媒体计算机(视频和声音等多媒体的压缩,转换和回放都是要依赖于浮点运算的),也正好随着 80486 的流行,逐渐普及开来。

在处理器上融合浮点运算依然是困难的。即使到今天,很多低端的处理器,都不带有浮点处理器。 所以,号称能够上天入地的,被移植到很多低端设备比如手机上的 Linux 内核,必然是不能支持浮点运算的,因为这样会破坏内核的可移植性。我们都知道, 在内核模式下,为了保证内核操作的原子性,一般在内核从事关键任务的时候所有中断是要被屏蔽的,用通俗的话说就是内核在做事情的时候,其他任何人不得打 扰。 如果内核支持浮点运算,不管是硬件实现也好,软件模拟也罢, 如果允许在内核中进行像浮点计算这样复杂而耗时的操作,整个系统的性能和实时响应能力会急剧下降。  即使是在硬件上实现的浮点运算,也不是件容易的事情,会耗费CPU较多的时钟周期,比如 Pentium 上的浮点数除法,需要耗费 39 个时钟周期才行,在流水线设计的CPU中,这种占用多个时钟周期的浮点运算会让整个流水线暂停,让CPU的吞吐量下降。在现代 CPU 设计中,工程师们发明了超标量,乱序执行,SIMD 等多种方式来克服流水线被浮点运算这种长周期指令堵塞的问题,这都是后话了。

正因为对于计算机来说,浮点运算是一个挑战性的操作,但又是做科学计算所需要的基本操作,所以浮点计算能力就成了计算机能力的一个测试标准。我们常常听说有一个世界上前 500 台最快的超级计算机列表,这里所谓的“快”的衡量标准,就是以每秒钟进行多少次浮点计算(FLOPS) 为准。按照 Top500.org, 即评选世界上前 500 台超级计算机的机构 2009年6月的数据,世界上最快的计算机,部署在美国能源部位于新墨西哥的洛斯阿拉莫斯国家实验室 (Los Alamos National Laboratory),当年造出第一颗原子弹的实验室。这台超级计算机,浮点计算速度的峰值高达 1456 TFlops,主要用来模拟核试验。因为美国的所有核弹头,海军核动力航母中的反应堆以及核试验,都由能源部国家核安全署(NNSA) 管理,所以能源部一直在投资用超级计算机进行核试验。 在 1996 年美国宣布不再进行大规模的物理核试验后的这么多年,美国能源部一直用超级计算机来做核试验,所以在 Top500 列表中,美国能源部拥有最多数量的超级计算机。

3. 数组下标寻址之障

言归正传,我们刚才说了在早期计算机发展史上,浮点计算的困难。除了浮点计算,还有一件事情特别困难,叫做数组下标寻址。用现代通俗的话说,就是当年的计算机,不直接支持 A[3] 这样的数组索引操作,即使这个操作从逻辑上说很简单:把数组 A 的地址加上 3,就得到了 A[3] 的地址,然后去访问这个地址。

这个困难在今天的程序员看来是不可思议的。 为什么这么简单的数组下标寻址机制最一开始的计算机没有支持呢? 原来,当年的计算机内存很小,只有一千到两千的存储空间,所以,描述地址只需要几位二/十进制数(BCD)。从而,在每条指令后面直接加一个物理地址是可行且高效的寻址方式。这种寻址方式,叫做直接寻址,当时所有的机器,都只支持直接寻址,因为在机器码中直接指出操作数的准确地址是最简单直接的方法,计算机不需要任何复杂的地址解码电路。但坏处是,这个设计太不灵活了,比如说 A[3] 这个例子,就没法用直接寻址来表示。

一般情况下,如果知道数组A, 对于 A[3] 这样的例子,用直接寻址问题去模拟间接寻址的困难还不是很大,只要程序员事先记住数组 A 的地址然后手工加上 3 就行了 (A也是程序员分配的,因为当时没有操作系统,所以程序员手工管理内存的一切)。可是,也有一些情况这样直接寻址是不行的。比如说,当时计算机已经能支持跳转和判断指令了,也就是说,可以写循环语句了。我们可以很容易看到, 以 i 为循环变量的循环体内,对 A[i] 的访问是不能写成一个静态的直接寻址的,因为 i 一直在变化,所以不可能事先一劳永逸的定好 A[i] 的所在位置,然后静态写在程序中。

这样,即使写一个简单的 10×10 矩阵的乘法,程序员就不得不死写 10的三次方即1000 行地址访问,而没办法用几行循环代替。当时的一些聪明人,也想了一些方法去克服这个问题,比如说,他们先取出 A 的地址,然后做一次加法,把结果,也就是当时 A[i] 的地址,注射到一个读内存的 LOAD 指令后面。然后执行那条 LOAD 指令。比如我要读 A[i],我先看,A的地址是 600,再看看 i 是3, 就加上 i,变成603,然后,把后面的指令改成 LOAD 603, 这样,就可以读到 A[i]。这个小技巧之所以可行,要感谢冯诺依曼爷爷的体系设计。在冯诺依曼计算机中,数据和程序是混在一起不加区分的,所以程序员可以随时像修改数据一样修改将要运行的下一条程序指令。就这样,靠着这个小技巧, 好歹程序员再也不要用1000行代码表示一个矩阵乘法了。

4. SpeedCoding 的出现

计算机本来就是用来做数学计算的,可是科学计算里面最最基本的两个要素–浮点计算和数组下标访问,在当时的计算机上都缺少支持。这种需求和实际的巨大落差,必然会召唤出一个中间层来消弭这种落差。 其实计算机科学的一般规律就是这样:当 A 和 C 相差巨大的时候,我们就引入一个中间层 B,用 B 来弥合 A 和 C 之间的不兼容。 当年的这个中间层,就叫做 SpeedCoding,由 IBM 的工程师 John Backus 开发。

SpeedCoding,顾名思义,就是让程序员编程更快。它其实是一个简单,运行在 IBM 701 计算机上的解释器。它允许程序员直接写浮点计算和下标寻址的指令,并且在底层把这些 “伪指令” 翻译成对应的机器码,用软件模拟浮点计算,自动修改地址等等。这样,程序员就可以从没完没了的手工实现浮点运算和下标寻址实现中解放出来,快速的编程。这个 SpeedCoding,这可以算得上是 FORTRAN 的种子了。

虽然这个解释器超级慢,程序员用这个解释器也用得很爽,也不感到它非常慢。 这是因为当年计算机浮点计算都绕不过软件模拟,即使最好的程序员用机器码而不用这个解释器,写出来的程序,也不比这个解释器下运行快多少。另一个更加重要的原因是,这个解释器极大的减少了程序员 debug 和 code 的时间。随着计算机速度的提高,当年一个程序耗费的计算成本和程序员编程耗费的人力成本基本上已经持平了,所以,相比较于写更加底层的机器码,用了 SpeedCoding 的程序员的程序虽然慢点,但人力成本瞬间降成 0,总体下来,用 SpeedCoding 比起不用来,总体成本还要低不少。

好景不长,因为客户一直的要求和电子工业的发展,IBM 在 1954 年,终于发布了划时代的 704 计算机,很多经典的语言和程序,都首次在 704 上完成了。比如之前我们在本系列的D篇中提到的 Steve Russell 的 LISP 解释器,就是在 704 上完成的。 704 计算机一下子支持了浮点计算和间接下标寻址。 这下用 SpeedCoding 的人没优势了,因为机器码支持浮点和下标寻址之后,写机器码比写 SpeedCoding 复杂不了多少,但是速度快了很多倍,因为 SpeedCoding 解释器太慢了,以前因为浮点和解释器一样慢,所以大家不在意它慢,现在浮点和寻址快了,就剩下解释器慢,写机器码的反而占了上风,程序员也就不用 SpeedCoding 了。

5. FORTRAN 创世纪

在 704 出来之前,做 SpeedCoding 的 John Backus 就认识到,要想让大家用他的 SpeedCoding, 或者说,想要从软件工具上入手,减少程序的开发成本,只有两个方法: 1. 程序员可以方便的写数学公式  2. 这个系统最后能够解析/生成足够的快的程序。他认为,只有达到了这两点,程序员才会乐意使用高级的像 SpeedCoding 这样的工具,而不是随着硬件的发展在机器码和 SpeedCoding 这样的工具之间跳来跳去。他本人通过实现 SpeedCoding, 也认识到如果有一个比机器码高级的语言, 生产效率会高很多倍。那么,现在唯一的问题就是实现它,当然,这就不是一个小项目了,就需要 IBM 来支持他的开发了。 所以,在 1953年,他把他的想法写成了一个文档,送给了 IBM 的经理。项目在 1954 年, 704 发布的当年,终于启动。John Backus 领导的设计一个能达到上面两点的编程系统的项目的成果,就是日后的 FORTRAN。

和现在大多数编程语言不一样,FORTRAN 语言的设计的主要问题不是语法和功能,而是编译器怎么写才能高性能。John Backus 日后回忆说,当时谁也没把精力放在语言细节上,语言设计很潦草的就完成了(所以其后正式发布后又经过了N多修订),他们所有的功夫都是花在怎么写一个高性能的编译器上。这个高性能的编译器很难写,到 1957 年才写好,总共花了 IBM 216 个人月。等到 FORTRAN 一推出,不到一年的时间,在 IBM 总共售出的 60 台 704上,就部署了超过一半。现在没啥编程语言能够这么牛的攻城掠地了 :)

6. 结语

放到历史的上下文中看,FORTRAN 的出现是很自然的。一方面,复杂的数学运算使得一个能够表述数学计算的高级语言成为必须,计算机的发展也为这个需求提供的硬件条件;另一方面,随着计算机的发展,程序员的时间成本一直不变,但是计算的成本一直在降低,用高级语言和用机器码在性能上的些许差异变得可以忽略。这样的历史现实,必然会召唤出以少量的增加计算机工作量为代价,但能大幅度降低程序员时间成本的新的工具和设计。这种新的工具,新的设计,又对程序设计产生革命性的影响。在整个编程发展的历史上,FORTRAN 和其他高级语言的出现可以说是第一波的革命;而后, UNIX和C语言的兴盛,使得系统编程的效率得到革命性提升,可以算是第二波革命;而面向对象方法,使得复杂的 GUI 等系统的编程效率得到提升,应该算得上是第三波革命。到如今,现在各种各样的方法论就更加多了,且看以后回看,哪种方法和工具能够大浪淘沙留下来。

(把八卦和文章重新分了一下类, 准备围绕自己的心得写个编程珠玑番外篇, 这个算第三篇吧)

<代码大全> (Code Complete) 是一本很好的书. 我建议像我这样写的程序总行数不超过50万的程序员应该买一本放在案头 (当然<代码大全> 不如 <Software Tools> 这本书好, 这个我以后有机会写文章细谈) 如果你天天编程序, 我建议你买一本. 如果你已经有了<代码大全>, 我诚心建议你赶快翻开此书, 撕去第26章. 因为代码大全的其他章节可以让你成为优秀的程序员, 唯独第26章, 读了之后立即从优秀程序员变成最差的程序员.

为啥? 因为第26章讲的, 都是怎么调节代码使得代码跑得更加快的技巧, 而这些技巧, 几乎都是让一个好程序变成差程序的技巧, 是教你不管三七二十一先对程序局部优化的技巧. 而局部优化是让程序变得糟糕的最主要的一个原因. 用高爷爷的话说, 提前优化是万恶之源 (Premature optimization is the root of all evil). 这些技巧, 就是带你去万恶之源的捷径.

代码优化究竟是什么洪水猛兽, 又究竟有多少伟大的程序员因为代码优化声名扫地, 请看本期关于代码优化的八卦.

话说当年在贝尔实验室. 一群工程师围着一个巨慢无比的小型机发呆. 为啥呢, 因为他们觉得这个机器太慢了. 什么超频, 液氮等技术都用了, 这个小型机还是比不上实验室新买的一台桌上计算机. 这些家伙很不爽, 于是准备去优化这个机器上的操作系统. 他们也不管三七二十一, 就去看究竟那个进程占用CPU时间最长, 然后就集中优化这个进程. 他们希望这样把每个程序都优化到特别高效, 机器就相对快了. 于是, 他们终于捕捉到一个平时居然占50% CPU 的进程, 而且这个进程只有大约20K的代码. 他们高兴死了, 立即挽起袖子敲键盘, 愣是把一个20K的C语言变成了快5倍的汇编. 这时候他们把此进程放到机器上这么一实验, 发现居然整体效率没变化. 百思不得其解的情况下他们去请教其他牛人. 那个牛人就说了一句话: 你们优化的进程, 叫做 System Idle.

所以说. 优化这东西, 一定要有一个全局的思路, 否则就是纯粹的无用功, 有时候还是负功. 在<编程珠玑 II> 第一章, Jon Bentley 就着重提醒了代码 profiling 的重要性. 说到 profiling 这个词, 就不能不再次提到万众敬仰的高爷爷. 高爷爷在1970年的暑假, 通过捡Stanford 大学机房扔出来的垃圾(其实是含有程序的磁带), 写出了一篇震古烁今的论文 “An empirical study of FORTRAN programs” (FORTRAN 程序的实证分析). 除了抱怨写程序的人不看他的 TAoCP 之外(因为一个程序用了被高爷爷定性为史上最差的随机数发生器算法, 有兴趣的可阅读 TAoCP vol2), 这篇论文主要说了三个划时代的东西:

1. 对程序进行 profile 是每个编程系统的居家旅行必备.

2. 在没 IO 操作的情况下, 一个程序中 4% 的代码占用了超过50% 的运行时间.

3. 97% 的情况下对程序进行提前优化是万恶之源.

这三个道理, 用大白话说, 就是: 1 程序都存在热点, 有优化的空间. 2. 但是97%的情况下程序员优化的都是错的地方, 反而把程序优化糟了. 3. 想要做优化, 第一步就要先知道程序在什么地方耗时间而不是靠猜.

说到热点, 顺带拐八卦一下Java的速度. Java 1.5 的虚拟机的关键技术, 就是叫做 Hotspot (热点). 传统上, 大家都认为 Java 比C 要慢. 其实不然. Jython 的作者 Jim Hugunin 就曾经说过, 其实两者差别不大 (http://hugunin.net/story_of_jython.html). 也有一些其他的测评说, Java 比 C 要快. 原因就在于, Java 虚拟机能够找到热点, 对热点专门做优化. 而C程序编译好了, 即使有热点, 也只能靠CPU去优化了. Java 的优化比 CPU 要深且更全局.

言归正传. 关于 FORTRAN 的 profile 的传统被继承了下来, 基本上现在任何的过程式主流编程语言都支持 profiling 工具. 关于 profile 怎么做的问题, 等我有空了好好写文章介绍. (因为我发现, 除了编程珠玑, 没有一本书提到过).

做程序优化的八卦就太多了, 说一个Beautiful Code 上的吧. 话说世界上做线性代数的库叫做BLAS, 基本上是工业标准. 因为线性代数运算太重要了, 所以各大处理器厂商都有 BLAS 的实现. Intel 的叫 MKL, AMD 的叫 ACML. 矩阵乘法实现的好坏, 直接决定了处理器的性能测试的分数(因为现代测处理器的速度的程序, 比如LAPACK指数, 基本上都是用 BLAS 里的矩阵乘法做基准). 去年 nVIDIA 高调宣传自己的 CUDA 系统比CPU厂商快10倍到100倍, 借此打开了GPU计算的大门(令人发指的达到500GFlops, Intel 最新的只有50GFlops). 其中 CUDA 可以理解为是 BLAS 在 nVIDIA 平台上的实现. 自从nVidia 推出 CUDA 以后, 俨然不把 intel 这些厂商放在眼里, 心想, 小样, 你们还是做通用处理器吧, 浮点乘法这些高级的东西, 还是放在显卡上比较好. nVidia 和 IBM/SONY 的阴谋很不小呢. 要是浮点计算比 Intel 快这么一两个数量级, 以后世界上前五百名超级计算机就全部变成什么用光纤网连起来的 PS3 机群, nVidia 显卡机群之类的. 人家外行见了计算机科学家肯定要问: 你们搞高性能计算机究竟是搞计算还是打网络游戏啊??

还是言归正传(我怎么老走题?), 简要的说一下 BLAS 优化. 单处理器上对 BLAS 的优化主要体现在对 cache 的高效使用. 矩阵乘法中, 如果矩阵都是按照行存储, 则在A*B中, 对B的访问是按列的. 假设B一行有N个元素, 那么在存储器中, 两个同列不同行的元素所在的存储单元相差N. 因此, 对B的访问并不是局部化的. 因为访问不局部化, 所以每次乘法, 都需要从内存中调一个 cache 单元到CPU. 这个极大的降低了处理器的执行速度. 因此, 矩阵乘法的优化的核心, 在于局部化B的访问. 反过来, 如果矩阵按照列存储, 则要局部化对A的访问. 关于怎样局部化访问还能获得正确的乘法, Beautiful Code 一书的第14章有非常好的讲解, 我就不废话了. 总之, 矩阵乘法局部化的好坏, 取决于一个机器的 cache 的大小.

多处理器和向量处理器就又不一样了. 要想利用好, 就要把计算任务平均的独立的分到不同的处理器上. 所以, 在这里, 优化就变成了分解成若干的小问题. 因此, 分治算法成了主流. 具体矩阵乘法怎么分块, 大学数学都讲了, 我就不废话了. 事实上, 之所以 nVidia 能和 Intel 干, 就是因为显卡上令人发指的有 64 个计算核心, 而 Intel 最牛X的才 4个. 那为啥 Intel 自己不多做几个核心呢? 因为 Intel 自己把自己带沟里去了 — Intel 处理器太复杂支持的功能太多了, 一块硅片上根本放不下很多核心. 而 nVidia 一直就是专用处理器, 每个核心功能简单, 可以做到很小.

Beautiful Code 第14章就是讲了随着计算机体系结构的变化, BLAS 是怎么进化的. Tanenbaum 曾经说过, 随着一个科技的出现, 某个 idea 可能就销声匿迹了. 但是说不定下一波科技再来的时候, 这个 idea 又复活了. BLAS 从串行, 到向量, 再到串行(带cache的RISC), 再到向量(Cell), 就是一个绝好的例子. 对这个进化史感兴趣的读者不可错过这一章妙文.