1 /* Copyright (c) 2015 Gerald Knizia
2 *
3 * This file is part of the IboView program (see: http://www.iboview.org)
4 *
5 * IboView is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * IboView is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with bfint (LICENSE). If not, see http://www.gnu.org/licenses/
16 *
17 * Please see IboView documentation in README.txt for:
18 * -- A list of included external software and their licenses. The included
19 * external software's copyright is not touched by this agreement.
20 * -- Notes on re-distribution and contributions to/further development of
21 * the IboView software
22 */
23
24 #include "Ir.h"
25 #include "IrAmrr.h"
26 #include "CxMemoryStack.h"
27 #include "CxAlgebra.h"
28 #include <algorithm> // for std::min.
29 #include <stdexcept> // for std::runtime_error
30
31 namespace ir {
32
33 void EvalSlmX_Deriv0(double *IR_RP Out, double x, double y, double z, unsigned L);
34
35
36 // For a given 3x3 unitary rotation matrix R[3][3], construct the matrix T[lc,l'c'] which makes
37 // S[l,c](R r) = \sum_{c',l'} T[lc,l'c'] S[l',c'](r)
38 // for all 3-vectors 'r' (i.e., compute the transformation between solid harmonic components
39 // induced by rotations in real space).
40 //
41 // Notes:
42 // - The matrix elements T[lc,l'c'] for l != l' vanish.
43 // - T[lc,l'c'] is stored at pOut[iSlcX(l,c) + nSlmX(MaxL) * iSlcX(l',c')].
44 // nStride is set to nSlmX(MaxL); pOut is allocated on Mem by this routine.
45 // - Order of trafo not checked. If it doesn't work, use transpose(T) instead of T.
EvalSlcXRotationTrafo(double * & pOut,size_t & nStride,unsigned MaxL,double const * R,ct::FMemoryStack & Mem)46 void EvalSlcXRotationTrafo(double *&pOut, size_t &nStride, unsigned MaxL, double const *R, ct::FMemoryStack &Mem)
47 {
48 nStride = nSlmX(MaxL);
49 Mem.Alloc(pOut, nStride*nStride);
50 void
51 *pFreeMe = Mem.Alloc(0);
52 unsigned
53 nSl = nSlmX(MaxL),
54 nCa = nCartX(MaxL);
55 double
56 *pA, *pB;
57 Mem.Alloc(pA, nSl * nCa);
58 Mem.Alloc(pB, nSl * nCa);
59
60 // while this could probably be done by employing some funky solid harmonic
61 // transformation theorems, we here just set up an equation system such that
62 //
63 // a_i = T * b_i,
64 //
65 // where a_i = Slc(R r_i) and b_i = Slc(r_i) for a sufficient basis
66 // of vectors r_i (with Slc components in the rows).
67 //
68 // The equation system is then linearly solved via LAPACK.
69 //
70 // Such a set of vectors is given, for example, by the cartesian vectors
71 // r_i = (x_i,y_i,z_i)
72 // with all x_i + y_i + z_i <= MaxL (for all l at once)
73 // or x_i + y_i + z_i = MaxL (for individual l)
74
75 // setup the matrices B = (b_i) and A = (a_i).
76 unsigned
77 iComp = 0;
78 for (unsigned i = 0; i <= MaxL; ++i)
79 for (unsigned j = 0; j <= i; ++j)
80 for (unsigned k = 0; k <= j; ++k) {
81 // make trial vector x and it's rotated counterpart R x.
82 double
83 r[3] = {double(i), double(j), double(k)},
84 Rr[3] = {0., 0., 0.};
85 for (unsigned a = 0; a < 3; ++ a)
86 for (unsigned b = 0; b < 3; ++ b)
87 Rr[a] += R[a + 3*b] * r[b];
88 // evaluate SlcX(x) and SlcX(R x).
89 EvalSlmX_Deriv0(&pB[iComp * nSl], r[0], r[1], r[2], MaxL);
90 EvalSlmX_Deriv0(&pA[iComp * nSl], Rr[0], Rr[1], Rr[2], MaxL);
91 iComp += 1;
92 }
93 assert(iComp == nCa);
94
95 // solve for T such that A = T B.
96 // By transposing both sides we get this into standard form:
97 // A^T = B^T T,
98 // where we solve for T. Due to LAPACK restrictions we do the
99 // transposition explicitly here.
100 double
101 *pAt, *pBt, *pSig, *pWork;
102 Mem.Alloc(pAt, nCa * nSl);
103 Mem.Alloc(pBt, nCa * nSl);
104 Mem.Alloc(pSig, nCa);
105 for (unsigned iSl = 0; iSl < nSl; ++ iSl)
106 for (unsigned iCa = 0; iCa < nCa; ++ iCa) {
107 pAt[iCa + nCa*iSl] = pA[iSl + nSl*iCa];
108 pBt[iCa + nCa*iSl] = pB[iSl + nSl*iCa];
109 }
110
111 // FIXME: do at least the SVDs individually per l. Results should be the same,
112 // but this would make it faster.
113 FORTINT
114 nRank = -1,
115 info = 0,
116 // minimal work space required according to docs.
117 nWork = 3*std::min(nCa,nSl) + std::max(std::max(2*std::min(nCa, nSl), std::max(nCa, nSl)), nSl);
118 nWork += 2*nCa*nSl*nSl; // <- we don't really care...
119 Mem.Alloc(pWork, nWork);
120 DGELSS(nCa, nSl, nSl, pAt, nCa, pBt, nCa, pSig, 1e-10, &nRank, pWork, nWork, &info);
121 if (info != 0)
122 throw std::runtime_error("MakeRotatedSlmXTransform: dgelss failed.");
123 if (nRank != nSl)
124 throw std::runtime_error("MakeRotatedSlmXTransform: equation system solution went wrong.");
125
126 // copy solution to output matrix.
127 for (unsigned iSl = 0; iSl < nSl; ++ iSl)
128 for (unsigned jSl = 0; jSl < nSl; ++ jSl)
129 pOut[iSl + nStride * jSl] = pBt[iSl + nCa * jSl];
130
131
132 // for (unsigned l = 0; l <= MaxL; ++ l) {
133 // double
134 // *IR_RP pT = &pOut[what?].
135 // (use addressing as above.. just make svds separately)
136 // }
137 Mem.Free(pFreeMe);
138 }
139
140
141 } // namespace ir
142