2016年11月15日星期二

Windows上如何在进程间传递socket句柄

问题的背景是这样的:我需要写一个程序A,来launch另一个程序B。B是一个server程序,它需要监听一个TCP端口。这个端口是A来分配的,A需要确保这个端口没有被占用。

我最初的办法是这样:写一个函数,遍历现在的tcp table,来判断哪些端口已经在使用中。

std::vector<int> getFreePorts(int begin, int end)
{
    char errbuf[256];
    if (end <= begin) {
        THROW_EXCEPTION("illegal argument");
    }
    std::unique_ptr<MIB_TCPTABLE2, decltype(std::free)*> pTcpTable(nullptr, std::free);
    ULONG ulSize = sizeof(MIB_TCPTABLE);
    for (int iter = 0; iter != 2; ++iter) {
        pTcpTable.reset((MIB_TCPTABLE2*)malloc(ulSize));
        ULONG dwRetVal = GetTcpTable2(pTcpTable.get(), &ulSize, FALSE);
        if (dwRetVal == NO_ERROR) break;
        if (dwRetVal != ERROR_INSUFFICIENT_BUFFER) {
            THROW_EXCEPTION("GetTcpTable2 failed,error = %ul", dwRetVal);
        }
    }
    std::vector<uint8_t> avail(end - begin, 1);
    for (int i = 0; i != pTcpTable->dwNumEntries; ++i) {
        auto& table = pTcpTable->table[i];
        int off = table.dwLocalPort - begin;
        if (off >= 0 && off < avail.size()) {
            avail[off] = 0;
        }
    }
    std::vector<int> ret;
    for (int i = 0; i != (int)avail.size(); ++i) {
        if (avail[i]) ret.push_back(i + begin);
    }
    return ret;
}
但是很奇怪,子进程bind有时候还是会fail。于是我就开始我的第二个版本:父进程创建socket并bind,子进程去listen并且accept。但是这时候就冒出来一个问题:windows的handle是指针,而不是int。指针怎么跨进程传递……出乎意料的是,它直接强转成size_t来传递。

于是父进程这样创建socket:
 SOCKET ListenSocket = ::socket(AF_INET, SOCK_STREAM, IPPROTO_TCP);
 if (ListenSocket == INVALID_SOCKET) {
  THROW_EXCEPTION("socket function failed with error: %u\n", WSAGetLastError());   
 }
 // The socket address to be passed to bind
 sockaddr_in service;
 service.sin_family = AF_INET;
 service.sin_addr.s_addr = inet_addr("0.0.0.0");
 service.sin_port = htons(port);
 int iResult = bind(ListenSocket, (SOCKADDR *)&service, sizeof(service));
 if (iResult == SOCKET_ERROR) {
  int err = WSAGetLastError();
  closesocket(ListenSocket);
  THROW_EXCEPTION("bind failed with error %u\n", err);   
 }
 return ListenSocket;
把上面的代码放入一个for循环中,然后把得到的socket通过命令行参数传递给子进程。在调用CreateProcess函数时,bInheritHandles参数应该是TRUE。如果是用java的process类创建子进程,那么不用做什么特别设置,bInheritHandles一直是TRUE。在子进程启动之后,父进程应当关闭这个socket。

子进程这样使用:
#include <WinSock2.h>
#include <Windows.h>
#include <stdio.h>
int main(int argc,char* argv[]) {
 WSADATA wsaData;
 WSAStartup(MAKEWORD(2, 2), &wsaData);
 FILE* fd=fopen("server.log","w");
 SOCKET ListenSocket = _strtoi64(argv[1],nullptr,10);
 //----------------------
 // Listen for incoming connection requests 
 // on the created socket
 if (listen(ListenSocket, SOMAXCONN) == SOCKET_ERROR) {
  int err = WSAGetLastError();
  closesocket(ListenSocket);
  fprintf(fd,"listen function failed with error: %d\n", err);
  return -1;
 }
 while (true) {
  SOCKET AcceptSocket = accept(ListenSocket, NULL, NULL);
  if (AcceptSocket == INVALID_SOCKET) {
   fprintf(fd, "accept failed with error: %ld\n", WSAGetLastError());
   fflush(fd);
   break;
  }
  else {
   fprintf(fd, "Client connected.\n");
   fflush(fd);
  }
 }
 fprintf(fd, "shutdown");
 fflush(fd);
 closesocket(ListenSocket); 
 fclose(fd);
 WSACleanup();
 return 0;
}

比较有趣的是,如果通过netstat这样的命令查看,会显示是父进程在监听端口并接受连接。这是一个绕过防火墙的很好的办法。Windows用户喜欢把一些程序列为可信任的,这些程序(比如IE、QQ)不受防火墙规则的限制。那么你只需要往这样的程序注入一小段代码,就可以让任意的程序悄悄绕过防火墙。这比直接修改防火墙规则要隐蔽多了。

上述方法适用于父进程和子进程存在父子关系的情况。如果不是父子关系,那么就得使用WSADuplicateSocket函数。source进程要先获知target进程的pid,然后通过WSADuplicateSocket函数得到一个WSAPROTOCOL_INFO结构体,然后自己想办法把WSAPROTOCOL_INFO这个结构体传递给target进程。

2016年11月14日星期一

进程内的生产者/消费者队列,一点思考和总结

考虑这样一个实际场景:
有一个文本文件,需要根据每一行的内容计算出一个score来。将score输出到一个新的文件。假设计算score是很耗费cpu的,请尽可能用多线程并行化的方式来加速。

这时我们会想到map-reduce/work-stealing/bounded queue等等。

拿bounded queue来说,最基本的版本是一个mutex两个条件变量,高级的版本是两个mutex两个条件变量。

接下来要解决的问题是:当生产者读到文件末尾时,如何把这件事情告诉消费者? 假设只有一个生产者,有多个消费者。一种简单的做法是:如果消费者的数量固定,已知,假设有n个,那么就往队列里扔n个特殊类型的消息以标记EOF。当消费者读到EOF则退出。这个思路也可扩展到多个生产者、多个消费者的情况。这是一种纯message-passing的方案。

但是如何处理突发出错?假设生产者读文件的时候突然读到一条错误的记录,它想要立刻通知所有的消费者结束,该怎么办?如果我们的队列是支持优先级的,那么把control message设置成最高优先级,这个问题也可以解决。但如果是反过来,消费者报错,希望所有的生产者和消费者都立刻终止,该怎么办?如果我们有一个双向的queue,这个问题也可以解决。但是往往最简单的方案还是加入一个coordinator,让它具有一个isHealthy()这样的方法。此时,所有条件变量的wait方法都需要改成timed wait。

再看另一个与此很相关的问题:假设我们有n个task,要交给n个线程同时执行。然后主线程需要等待并搜集这n个结果。这个看似很简单,可以用count down latch实现一个栅栏,也可以创建n个feature 然后挨个wait。但是假如要考虑出错的情况,即,子线程在执行的时候可能会遇错提前退出,此时希望其它的线程也尽快停止。

这时候最简单的办法是:让count down latch有一个countDownToZero这样的方法。这样主线程就收到错误直接停止等待了。还有一种办法就是把task result及task id都放入一个unbounded queue中,然后主线程一直wait在这个queue上。

但无论如何,为了能安全且及时的通知其它线程退出,无论什么样的场景,代码中不能有无限时间的wait。所有的wait必须得带有一个时间间隔。相比而下java就简单些,因为java可以interrupt线程。C/C++中虽然有pthread_cancel这样的方法,但是终归还是太危险。

事情还没完,假如IO是异步的,这问题又复杂许多。举个例子,假设读文件是先发起一个http请求,然后等待callback处理IO完成或异常。那么怎么处理出错呢?我们可以在每个callback的开始加上isHealthy(),此外,还必须知道有没有callback正在执行。如果有,主线程就得继续等。如果没有,必须保证不会有新的callback被触发,然后才能退出。

2016年10月27日星期四

如何在VC++中强制引用静态库里的所有全局变量

这个坑啊…… 真是历史悠久名声昭著了。

假如你的exe用到了一个静态库,而这个静态库有一些全局变量。那么这些全局变量未必会被引用到这个exe当中。

举个例子:
Foo.h:
===================
class Foo
{
public:
Foo();
};
===================

Foo.cpp:
===================
#include "Foo.h"
#include <stdio.h>


Foo::Foo()
{
printf("hello");
}

static Foo foo;
===================
假如Foo.cpp是被编译成静态库,那么链接的时候可能没有foo这个变量,从而也不会有"hello"输出。

有些人写代码喜欢用全局变量来实现singleton,并且期望这些singleton的构造函数会在main函数之前被调用。呵呵…… 到了vc这里就不灵了。

可是,如果万一他的代码就是这样,你又没法改,办法还是有的。

1. 要以项目引用的方式来引用静态库,而不是把.lib文件作为input给linker
2. 把“link library dependencies” 和 "Use link library dependency inputs” 这两个选项设置成true。

亲测有效。VC++ 2015。


2016年10月22日星期六

如何正确的关闭已打开的文件

如果一个文件是以只读的方式打开的,那么忘记关闭的后果是:句柄泄露。通常来说这不是太大的问题。
如果一个文件是以可写的方式打开的,那么忘记关闭的后果是:除了句柄泄露以外,写进去的东西可能会丢失。并且,如果关闭了,但是忘记检查close(或fclose)函数的返回值,也会有丢数据的风险。

要注意点什么?

写文件的时候,要手动关,不要依赖于析构函数。用RAII来管理IO资源是非常愚蠢的。如果fclose出现在析构函数里,那么通常就是一个错误。因为析构函数不能抛异常,如果fclose返回非0的值,你很难把这个错误报告出去。出错不可怕,最可怕的是出错了但是你却不知道。比如,你的程序每天自动生成一个config文件然后自动推送到线上每台机器。结果有一天,生出来的config文件缺了一块,但是你却不知道,而这个config却被推送出去了。

什么时候fclose会返回非0值?

比如,硬盘满了。再比如,有些机器的硬盘会有一些临时故障,会导致写入偶尔失败。而这种临时故障不一定会产生报警。其实我们希望一旦有这样的问题发生,这些机器能赶紧被应用程序发现,然后从系统中剔除出去。

正确的该怎么写?

FILE* fd = fopen("xxx.txt","w");
if(!fd) return;
try{
  //do normal io stuffs
   ... 
}catch(std::exception& ex){
   if(fclose(fd)){
      //....?
   }
}
if(fclose(fd)) throw std::runtime_error("err");

问题来了,如果我们想把fclose失败这件事情报告出去,就得在处理异常的时候又抛出新的异常。该怎么办?遗憾的是,C++中对于这样的问题,并没有标准的做法。

Java是怎么处理这样的问题?

Java7为此特地推出了Suppressed Exceptions.
比如下面的代码:
package testjava;

import java.io.Closeable;
import java.io.IOException;

public class T1 {

 static class A implements Closeable{

  @Override
  public void close() throws IOException {
   throw new RuntimeException("Err");   
  }
  
 }
 
 static void foo() throws Exception{
  try(A a=new A()){
   throw new RuntimeException("er");
  }
 }
 
 public static void main(String[] args) throws Exception {
  try{
   foo();
  }catch(Exception ex){
   ex.printStackTrace();
  }
 }

}
会产生下面这样的输出:

java.lang.RuntimeException: er
at testjava.T1.foo(T1.java:19)
at testjava.T1.main(T1.java:25)
Suppressed: java.lang.RuntimeException: Err
at testjava.T1$A.close(T1.java:12)
at testjava.T1.foo(T1.java:20)
... 1 more

虽然代码看起来更简单了,但是这其实是把问题隐藏的更深了。谁在处理异常的时候会检查它内部有没有另一个异常呢?所以你其实没有机会把硬盘错误报告出去。

最终最重要的一句话:
“Don't trust your callers, nor your callees, always safeguard your code”

2016年10月6日星期四

人民币与IMF特别提款权(SDR)

最近关于人民币被纳入IMF的特别提款权的讨论很热,我认同这个观点:“它的象征意义大于实际意义。”

SDR是布雷顿森林体系的产物,在2009年以前从未发挥过实际作用。在布雷顿森林体系时代,各国央行为了维持本国货币与美元之间的固定利率,不得不储备大量美元。于是SDR就被设计出来,作为外汇储备的一种。可惜它还尚未发挥实际作用,固定汇率制度就倒台了,于是SDR也就被尘封起来。然而若干年后,中国崛起,中国央行长久以来选择了人民币和美元之间的固定汇率制度。如下图:

1美元兑人民币
为了干预汇率以及强制结汇,中国央行不得不储备大量的美元,代价是每年得为此付出高昂的成本。2008年之后,美元大幅贬值,但中国央行继续增加美元外汇储备大量购买美国国债,于是遭到很多质疑。于是央行行长周小川就瞄准了SDR。2009年3月,周小川在中国人民银行网站上发文《关于改革国际货币体系的思考》,文中建议“创造一种与主权国家脱钩、并能保持币值长期稳定的国际储备货币,从而避免主权信用货币作为储备货币的内在缺陷”。这种货币指的就是SDR。它希望在外汇储备资产中,增加SDR的比重,减少美元,以降低美国对中国的货币经济影响。

但SDR不是你想买就能买的,它的发行量非常小。当时SDR在全球的发行量只有20亿美元,中国拥有4%,大约8千万,和其上万亿的外汇储备相比不值得一提。为了提高SDR在外汇储备中的占比,中国做了两件事情:

  1. 推动IMF在2009年对SDR做了一次大规模增发,使其发行量一次性扩大10倍左右。
  2. 要求增缴份额。使其在IMF份额从3.996%升至 6.394%,投票权随之增加。其SDR allocation在将来也会随之同比扩大。
  3. 要求把RMB纳入SDR的货币篮子当中
第一条,不需要美国国会批准。而对于后两点,美国是极力反对的。所以这件事情从2010年被提出,一直拖到了2015年底,才终于在美国国会批准。可以说,不是IMF批准人民币“入篮”SDR,而是美国国会批准了此事。因为美国对IMF具有绝对的控制权,它在IMF中有权否决一切提议。

截至2016年6月30日,中国央行的外汇储备资产中包括:3万亿的外汇储备、800亿的货币黄金和100亿的SDR。可见SDR的地位可忽略不计,更何况SDR的流动性远远不能和美国国债相比。

有观点说,SDR作为一种外汇储备资产,可被列入IMF 100多个成员国的资产负债表上。而RMB作为SDR的一部分,储备SDR就相当于储备RMB,所以RMB现在是100多个国家的外汇储备货币。呵呵…… 这么看象征意义是蛮大的。

另外一件很具象征意义的事情是:世界银行在中国银行间债券市场发行首期特别提款权(SDR)计价债券,发行规模为5亿SDR(约合46亿RMB),期限为3年,结算货币为人民币。这是世界银行发行的第一个以SDR计价的债券。票面利率仅为0.49%。按理说利率这么低应该没人买,结果被一抢而空。据推测原因主要是:所有人都预期人民币会继续贬值。今年上半年RMB相当于SDR已经贬值3%,如果照此速度,这个债券的票面利率将高达20%。呵呵……当然实际不会这么狠啦。


2016年9月26日星期一

How to removes duplicate values from an array: A counterintuitive story of Big O

Everyone knows O(N) is faster than O(NlogN). However, sometimes, it isn't.
Please look at this example:

template <typename T>
size_t count_unique_elements_by_map(T* begin, T* end) {
  std::unordered_set<T> s(end - begin);
  s.insert(begin, end);
  return s.size();
}

template <typename T>
size_t count_unique_elements_by_sort(T* begin, T* end) {
  if (begin >= end) return 0;
  std::sort(begin, end);
  auto newEnd = std::unique(begin,end);
  return newEnd - begin;
}

Method 2 is faster than method 1 if duplicates are rare.

问题的背景是这样:我有一些ID,我想把这些ID对应的value查出来。这些value保存在一个分布式系统中(如memcached)。为了减少服务器的查询压力,我把这些ID先去重,然后计算每个ID在哪个服务器上,把ID按服务器分组,然后发送查询。最初我用了一个hash map做这件事情,然后我发现去重的操作所耗费的CPU远远超过其它代码(比如网络IO)。

然后我发现sort比hash更快,在我的场景下,快10倍以上。但是故事并未到此结束,当我收到服务器的reply之后,我还是得把这些ID以及value插入到一个hash map中。否则我就只能用binary search读取它们。所以,我还是得构建这个hash map,逃不掉的。

再后来我发现,如果把STL的hash map换成我自己写的一个基于open address方式的hash map,那么速度也能提升1倍左右。

Please tell me if you have any better solutions. Thanks.

update: 我刚想到,如果input id也是有序的,那么不用构建hash map。在有序数组中查找有序数组,算法复杂度还是O(N)。但是实际效果如何,有待验证。

Torch的tensor

Torch中的tensor是一个n维数组,它是torch最基础最重要的组件。
下面这段代码演示如何通过下标来访问tensor。这段代码首先,初始化一个长度为10的1维数组,将其内容填充成0,1,2...9,然后遍历这个tensor,将其内容打印出来。
#include "TH/TH.h"

#define Real Double
#define real double

int main() {
 THTensor* tensor1 = THTensor_(newWithSize1d)(10);
 for (int64_t i = 0; i != 10; ++i) {
  THTensor_(set1d)(tensor1, i, (real)i);  
 }
 for (int i = 0; i != 10; ++i) {
  printf("%g ", THTensor_(get1d)(tensor1, i));
 }
 printf("\n");
 THTensor_(free)(tensor1);
 return 0;
}
输出应该是:
0 1 2 3 4 5 6 7 8 9
使用下标来访问tensor有两种方式,上面这种是通过set/get函数,比较慢。更快的方式是这样。
#include "TH/TH.h"

#define Real Double
#define real double

int main() {
  THTensor* tensor1 = THTensor_(newWithSize1d)(10);
  int i = 0;
  TH_TENSOR_APPLY(real, tensor1, *tensor1_data = (real)i; ++i;);
  for (int i = 0; i != 10; ++i) {
    printf("%g ", THTensor_(get1d)(tensor1, i));
  }
  printf("\n");
  THTensor_(free)(tensor1);
  return 0;
}
TH_TENSOR_APPLY是一个宏,在TH/THTensorApply.h中声明,它的code很复杂,但是又很需要读一读。本文之后的内容主要介绍这个宏的原理。

首先看下tensor的内部构造。
struct THTensor
{
    long *size;
    long *stride;
    int nDimension;
    
    THStorage *storage;
    long storageOffset;
    int refcount;
    char flag;
}
tensor的实际内容存在THStorage中。拿数据库来做类比,如果把THStorage看成是物理表,那么tensor就是视图。下面是从tensor获取storage指针,然后从storage通过THStorage_(data)函数获得裸指针,然后遍历这个tensor。如下面这段代码所示。
#include "TH/TH.h"

#define Real Double
#define real double

int main() {
 THTensor* tensor1 = THTensor_(newWithSize1d)(10);
 THStorage* s = THTensor_(storage)(tensor1);
 double* data = THStorage_(data)(s);
 int64_t size = THStorage_(size)(s);
 for (int64_t i = 0; i != size; ++i) {
  data[i] = (double)i;
 }
 THTensor_(free)(tensor1);
 return 0;
}
所以当需要实现向量的加减乘除、点乘等操作时,出于效率考虑,肯定会采用后一种方法,而不是采用set1d/get1d。于是问题来了。给定下标,如何计算出offset ?
以一个4维的tensor为例子
tensor[d][c][b][a] 的offset等于  a*stride[0] + b*stride[1] + c *stride[2] + d*stride[3]]
如下面的代码所示:
#include "TH/TH.h"
#include <assert.h>

#define Real Double
#define real double

int main() {
 THTensor* ts = THTensor_(newWithSize4d)(8, 4, 6, 7);
 long stride1 = THTensor_(stride)(ts, 0);
 long stride2 = THTensor_(stride)(ts, 1);
 long stride3 = THTensor_(stride)(ts, 2);
 long stride4 = THTensor_(stride)(ts, 3);
 printf("strides: %ld %ld %ld %ld\n", stride1, stride2, stride3, stride4);
 double* data = THTensor_(data)(ts);
 int a = 2;
 int b = 3;
 int c = 5;
 int d = 2;
 THTensor_(set4d)(ts, a, b, c, d, 12);
 assert(data[a*stride1 + b*stride2 + c *stride3 + d*stride4] == 12);
 THTensor_(free)(ts);
 return 0;
}
输出:
strides: 168 42 7 1
THTensor_(data) 函数返回的并不是struct THTensor中storage的data指针,而是data+storageOffset。

real *THTensor_(data)(const THTensor *self)
{
  if(self->storage)
    return (self->storage->data+self->storageOffset);
  else
    return NULL;
}

大多数时候我们无需关注这个细节。

对于刚new出来的tensor,torch可以保证tensor的最后一维总是连续的。所以strides数组的最后一个元素一定是1。它的strides满足:
stride[d] = isLastDimension(d) ? 1 : size[d+1] * stride[d+1];
但情况并非总是如此。比如一个2维的tensor经过转置之后,它依然和之前的tensor共享同一个storage。比如下面这段代码:
#include "TH/TH.h"
#include <assert.h>
#include <stdexcept>

#define Real Double
#define real double

static void printMatrix(THTensor* tensor) {
 int nDim = THTensor_(nDimension)(tensor);
 if (nDim != 2)
  throw std::runtime_error("dimension error");
 long size0 = THTensor_(size)(tensor, 0);
 long size1 = THTensor_(size)(tensor, 1);
 for (int i = 0; i != size0; ++i) {
  for (int j = 0; j != size1; ++j) {
   printf("%g ", THTensor_(get2d)(tensor, i, j));
  }
  printf("\n");
 }
}

int main() {
 THTensor* tensor1 = THTensor_(newWithSize2d)(3,4);
 THStorage* s = THTensor_(storage)(tensor1);
 long stride1 = THTensor_(stride)(tensor1, 0);
 long stride2 = THTensor_(stride)(tensor1, 1);
 printf("Before transpose:\n");
 printf("strides: %ld %ld\n", stride1, stride2);
 double* data = THStorage_(data)(s);
 int64_t size = THStorage_(size)(s);
 for (int64_t i = 0; i != size; ++i) {
  data[i] = (double)i;
 }
 printf("content:\n");
 printMatrix(tensor1); 
 printf("After transpose:\n");
 THTensor* tensor2 = THTensor_(newTranspose)(tensor1, 0, 1); 
 stride1 = THTensor_(stride)(tensor2, 0);
 stride2 = THTensor_(stride)(tensor2, 1); 
 printf("strides: %ld %ld\n", stride1, stride2);
 printf("content:\n");
 printMatrix(tensor1);
 THTensor_(free)(tensor2);
 THTensor_(free)(tensor1);
 return 0;
}
输出:
Before transpose:
strides: 4 1
content:
0 1 2 3
4 5 6 7
8 9 10 11
After transpose:
strides: 1 4
content:
0 1 2 3
4 5 6 7
8 9 10 11


这时候,如果想要遍历这个转置后的矩阵,就复杂的多了。对于一个n维数组,如果我们想遍历它,可以分配一个counter
size_t counter[dim]; //当前元素的下标
for(int i=0;i!=dim;++i) counter[0] = 0;

然后每遍历一个元素后,用下面的代码递增counter
bool isFinished = false;
for(int i=dim;i>=0;i--){
   counter[i]++;
   if(counter[i] != size[i]) break;
   if(i==0) {
     isFinished = true;
     break;
  }
   counter[i] = 0;
}
当isFinished=true的时候,终止所有循环。
这又变成了刚才的问题,如果我们对每个下标都用一系列的乘法和加法来计算offset,就会变得很低效。于是我们可以这样,用real* tensor1_data指向当前元素,然后在每次修改counter数组的时候,顺便也修正tensor1_data。
    for(i = tensor1_dim; i >= 0; i--) 
    { 
      tensor1_counter[i]++; 
      tensor1_data += tensor1->stride[i]; 

      if(tensor1_counter[i]  == tensor1->size[i]) 
      { 
        if(i == 0) 
        { 
          hasFinished = 1; 
          break; 
        } 
        else 
        { 
          tensor1_data -= tensor1_counter[i]*tensor1->stride[i]; 
          tensor1_counter[i] = 0; 
        } 
      } 
      else 
        break; 
    }
在此基础上可以再有另外一个优化,计算tensor的最低维部分有多少维是连续的,即满足 stride[d] = isLastDimension(d) ? 1 : size[d+1] * stride[d+1];
对于这些维度,我们可以用一个for循环一次遍历完。
完整的代码如下:(How to loop an n-dim array)
  real *tensor1_data = NULL; 
  long *tensor1_counter = NULL; 
  long tensor1_stride = 0, tensor1_size = 0, tensor1_dim = 0, i; 
  int hasFinished = 0; 

  if(tensor1->nDimension == 0) 
    hasFinished = 1; 
  else 
  { 
    tensor1_data = tensor1->storage->data+tensor1->storageOffset; 

    /* what is the first stride (ignore first dims=1)? */ 
    /* it will be used for the whole largest contiguous section */ 
    for(tensor1_dim = tensor1->nDimension-1; tensor1_dim >= 0; tensor1_dim--) 
    { 
      if(tensor1->size[tensor1_dim] != 1) 
        break; 
    } 
    tensor1_stride = (tensor1_dim == -1 ? 0 : tensor1->stride[tensor1_dim]); 

    /* what is the largest contiguous section? */ 
    tensor1_size = 1; 
    for(tensor1_dim = tensor1->nDimension-1; tensor1_dim >= 0; tensor1_dim--) 
    { 
      if(tensor1->size[tensor1_dim] != 1) 
      { 
        if(tensor1->stride[tensor1_dim] == tensor1_size) 
          tensor1_size *= tensor1->size[tensor1_dim]; 
        else 
          break; 
      } 
    } 

    /* counter over found dimensions */ 
    tensor1_counter = (long*)THAlloc(sizeof(long)*(tensor1_dim+1)); 
    for(i = 0; i <= tensor1_dim; i++) 
      tensor1_counter[i] = 0; 
  } 

  while(!hasFinished) 
  { 
    for(i = 0; i < tensor1_size; i++, tensor1_data += tensor1_stride) 
    { 
      //PUT USER CODE HERE 
    } 

    if(tensor1_dim == -1) 
       break; 
 
    tensor1_data -= i*tensor1_stride; 
    for(i = tensor1_dim; i >= 0; i--) 
    { 
      tensor1_counter[i]++; 
      tensor1_data += tensor1->stride[i]; 

      if(tensor1_counter[i]  == tensor1->size[i]) 
      { 
        if(i == 0) 
        { 
          hasFinished = 1; 
          break; 
        } 
        else 
        { 
          tensor1_data -= tensor1_counter[i]*tensor1->stride[i]; 
          tensor1_counter[i] = 0; 
        } 
      } 
      else 
        break; 
    } 
  } 
  THFree(tensor1_counter);