{"action":"created","issue":{"url":"https://api.github.com/repos/numpy/numpy/issues/619","labels_url":"https://api.github.com/repos/numpy/numpy/issues/619/labels{/name}","comments_url":"https://api.github.com/repos/numpy/numpy/issues/619/comments","events_url":"https://api.github.com/repos/numpy/numpy/issues/619/events","html_url":"https://github.com/numpy/numpy/issues/619","id":7718635,"number":619,"title":"BLAS matrix product (dot) never used for ndim > 2 (tensordot does not use BLAS) (Trac #2163)","user":{"login":"numpy-gitbot","id":1955252,"avatar_url":"https://avatars.githubusercontent.com/u/1955252?v=3","gravatar_id":"","url":"https://api.github.com/users/numpy-gitbot","html_url":"https://github.com/numpy-gitbot","followers_url":"https://api.github.com/users/numpy-gitbot/followers","following_url":"https://api.github.com/users/numpy-gitbot/following{/other_user}","gists_url":"https://api.github.com/users/numpy-gitbot/gists{/gist_id}","starred_url":"https://api.github.com/users/numpy-gitbot/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/numpy-gitbot/subscriptions","organizations_url":"https://api.github.com/users/numpy-gitbot/orgs","repos_url":"https://api.github.com/users/numpy-gitbot/repos","events_url":"https://api.github.com/users/numpy-gitbot/events{/privacy}","received_events_url":"https://api.github.com/users/numpy-gitbot/received_events","type":"User","site_admin":false},"labels":[{"url":"https://api.github.com/repos/numpy/numpy/labels/component:%20numpy.core","name":"component: numpy.core","color":"ff0000"},{"url":"https://api.github.com/repos/numpy/numpy/labels/Enhancement","name":"Enhancement","color":"e102d8"},{"url":"https://api.github.com/repos/numpy/numpy/labels/priority:%20normal","name":"priority: normal","color":"ff0000"}],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":5,"created_at":"2012-10-19T15:10:06Z","updated_at":"2015-08-04T03:20:14Z","closed_at":null,"body":"_Original ticket http://projects.scipy.org/numpy/ticket/2163 on 2012-06-15 by trac user thatistosay, assigned to unknown._\n\ndotblas_matrixproduct() contains the comment \"This function doesn't handle dimensions greater than 2\" and calls PyArray_MatrixProduct2() for these cases. This means BLAS is never used for calls to dot() with arguments of ndim>2!! In more detail...\n\nIf I want to contract a pair of tensor indices for ndim=3 such that \n\n\n\n    A = np.rand(d,D,D); B = np.rand(d,D,D)\n    AB[s,i,t,j] == sum(A[s,i,:], B[t,:,j])\n\n\n\nthen currently, although it can be done in a single line\n    \n\n\n    res = sp.dot(A,B)\n\n\n    \nit can often be done much faster with explicit (python!) loops\n\n\n\n    res = np.zeros((d,d,D,D))\n    for s in xrange(d):\n        for t in xrange(d):\n            np.dot(A[s], B[t], out=res[s,t])\n    res = np.rollaxis(res, 2, 1)\n\n\n\n..assuming dot() is using optimized BLAS for ndim=2, and the dimensions are large enough so that calling BLAS is worth it. \nIn general, reproducing the behaviour of dot() for ndim>2 is just a matter of calling GEMM in loops as above and then calling rollaxis() once.\n\nI therefore propose doing this within _blasdot.c as far as possible (to eliminate the use of python loops) so that ndim>2 dot(), and tensordot(), can benefit from BLAS.\n\nSome comparisons of the two methods above (attached script):\n\n\n    dtype=complex128\n    AB[s,i,t,j] = sum(A[s,i,:], B[t,:,j])\n    \n    A.shape = (16, 512, 512); B.shape = (16, 512, 512)\n    looping over 2D dot() vs. 3D dot(): 24% (about 4 times faster)\n    \n    A.shape = (20, 64, 64); B.shape = (20, 64, 64)\n    looping over 2D dot() vs. 3D dot(): 35%\n    \n    A.shape = (20, 48, 32); B.shape = (20, 32, 48)\n    looping over 2D dot() vs. 3D dot(): 45%\n    \n    A.shape = (32, 32, 16); B.shape = (32, 16, 32)\n    looping over 2D dot() vs. 3D dot(): 82%\n    \n    A.shape = (64, 10, 8); B.shape = (64, 8, 10)\n    looping over 2D dot() vs. 3D dot(): 158% (slow python loops..)\n\n\n(this was on a 4-core i7 system using ATLAS under heavy load)"},"comment":{"url":"https://api.github.com/repos/numpy/numpy/issues/comments/127463380","html_url":"https://github.com/numpy/numpy/issues/619#issuecomment-127463380","issue_url":"https://api.github.com/repos/numpy/numpy/issues/619","id":127463380,"user":{"login":"ewmoore","id":2314641,"avatar_url":"https://avatars.githubusercontent.com/u/2314641?v=3","gravatar_id":"","url":"https://api.github.com/users/ewmoore","html_url":"https://github.com/ewmoore","followers_url":"https://api.github.com/users/ewmoore/followers","following_url":"https://api.github.com/users/ewmoore/following{/other_user}","gists_url":"https://api.github.com/users/ewmoore/gists{/gist_id}","starred_url":"https://api.github.com/users/ewmoore/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ewmoore/subscriptions","organizations_url":"https://api.github.com/users/ewmoore/orgs","repos_url":"https://api.github.com/users/ewmoore/repos","events_url":"https://api.github.com/users/ewmoore/events{/privacy}","received_events_url":"https://api.github.com/users/ewmoore/received_events","type":"User","site_admin":false},"created_at":"2015-08-04T03:20:14Z","updated_at":"2015-08-04T03:20:14Z","body":"Given that numpy dot does an outer product type thing when given ndim > 2, and I'm not sure how the code he links to does ndim != 2 since i don't think the scipy he is depending handles it. I'd say the approach outlined here is much superior. \r\n\r\n"}}