5
$\begingroup$

I need to populate a matrix $A_{kl}$, where

$$ k = (m-1)J+n$$ $$ l = (p-1)J+q$$

And

$$m,p = 1, 2, ..., I$$ $$n,q = 1, 2, ..., J$$

Its components are (mnpq). For populate it, I'm using an expensive 4 Do loop

Do[
 Do[
  Do[
   Do[
    Print[m, n, p, q];
    k = (m - 1) nC + n;
    l = (p - 1) nC + q;
    If[k <= l, A[[k, l]] = cf[Nfunc, xi, yi, wix, wiy, m, n, p, q], 
     0];
    , {q, 1, J, 1}]
   , {p, 1, I, 1}]
  , {n, 1, J, 1}]
 , {m, 1, I, 1}]

Knowing that $A_{kl}$ for a $I=J=2$, its components are (mnpq)

$$ \begin{bmatrix} (1111) & (1112) & (1121) & (1122)\\ & (1212) & (1221) & (1222)\\ symm. & & (2121) & (2122)\\ & & & (2222)\\ \end{bmatrix} $$

Does anyone know a more efficient way to populate it? Maybe using a bult-in function?

UPDATE

cf = Compile[{{Nfunc, _Real, 2}, {xi, _Real, 1}, {yi, _Real, 
     1}, {wix, _Real, 1}, {wiy, _Real, 
     1}, {m, _Integer}, {n, _Integer}, {p, _Integer}, {q, _Integer}},
   Module[{sum},
    sum = 0.0;
    For[i = 1, i <= Length@xi, i++,
     For[j = 1, j <= Length@yi, j++,
      sum = 
        sum + (8 \[Pi]^2)/
           a^2 m p Cos[(m \[Pi] xi[[i]])/((1/
                2) a)] Sin[(n \[Pi] yi[[j]])/((1/
                2) b)] Cos[(p \[Pi] xi[[i]])/((1/
                2) a)] Sin[(q \[Pi] yi[[j]])/((1/2) b)]*wix[[i]]*
          wiy[[j]]*Nfunc[[i, j]];
      ]
     ];
    sum],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"];
$\endgroup$
6
  • $\begingroup$ At least you can just use only one Do. $\endgroup$ – Αλέξανδρος Ζεγγ May 5 '20 at 13:32
  • $\begingroup$ What is the nature of cf? Is it Listable in any of its arguments? I.e. would cf[nFunc, xi, yi, wix, wiy, m, n, p, {q1, q2, q3}] correctly give a list of three values? $\endgroup$ – Marius Ladegård Meyer May 5 '20 at 13:35
  • $\begingroup$ @ΑλέξανδροςΖεγγ, that is the exact point. I'm not getting create a algorithm to only one Do $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 13:53
  • $\begingroup$ @MariusLadegårdMeyer, m,n,p,q are integers. All others are lists. cf = Compile[{{Nfunc, _Real, 2}, {xi, _Real, 1}, {yi, _Real, 1}, {wix, _Real, 1}, {wiy, _Real, 1}, {m, _Integer}, {n, _Integer}, {p, _Integer}, {q, _Integer}}. $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 13:55
  • $\begingroup$ Yes, but you want to evaluate this function not for only one value of these integers, but for many different ones. If we put these many different ones in a list, and feed that list to cf, will it work? For instance, if I want to find Sin[x] for x = 1, 2, 3, ..., 100 I don't need to loop, I just do Sin[Range[100]]. I'm asking whether this applies to your function cf. $\endgroup$ – Marius Ladegård Meyer May 5 '20 at 13:58
5
$\begingroup$

Your $m-1, n-1$ and $p-1, q-1$ are the two digits of $k-1$ and $l-1$, respectively, in base $J$. This should be pretty fast:

dim = i*j;
A = ConstantArray[0, {dim, dim}];
Do[
   A[[k, l]] = A[[l, k]] = cf[Nfunc, xi, yi, wix, wiy, Quotient[k - 1, j] + 1, Mod[k - 1, j] + 1, Quotient[l - 1, j] + 1, Mod[l - 1, j]+ 1]
, {k, 1, dim}, {l, k, dim}
]

If $IJ$ is huge and your cf is very very fast, then you can probably shave off a little bit more time by making an outer Do over k where you compute the quotient and mod once for that k, followed by an inner Do over l.

$\endgroup$
4
  • $\begingroup$ I will update my question adding cf. Could you analyse it? Thank you for answering! $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 14:29
  • $\begingroup$ It doing i = 2; j = 2; dim = i*j; A = ConstantArray[0, {dim, dim}]; Do[A[[k, l]] = A[[l, k]] = {Quotient[k, j] + 1, Mod[k, j], Quotient[l, j] + 1, Mod[l, j]}, {k, 1, dim}, {l, k, dim}], we do not get the @bill s answer, as should be. $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 14:58
  • $\begingroup$ I have fixed the off-by-one issue. $\endgroup$ – Marius Ladegård Meyer May 5 '20 at 15:10
  • $\begingroup$ Perfect, @Marius. Thank you very much! $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 15:29
5
$\begingroup$

One way to create your data structure is to realize that your list of integers is closely related to all the Tuples taken in groups. For the case I=J=2, the permutations are of 1 and 2 taken 4 at a time:

mat = Partition[Tuples[{1, 2}, 4], 4];
mask = UpperTriangularize[ConstantArray[1, {4, 4}]];
outMat = mat mask + Transpose[mat mask] - mat IdentityMatrix[4]

outMat//MatrixForm

enter image description here

$\endgroup$
4
  • $\begingroup$ Thank you answering. But I did not understand how I will access each component of the matrix A using your approach. $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 14:28
  • 1
    $\begingroup$ Use Part ([[ ]]). For example, outMat[[2, 3]] will give you the 2,3 element of the matrix. $\endgroup$ – bill s May 5 '20 at 15:07
  • $\begingroup$ What a dumb! Sorry for this dumb question! Thank you very much! I accepted the @Marius's answer because it was more straight for me. But, without doubt, your approach is really great too! $\endgroup$ – Professor P. Cosmo Klunk May 5 '20 at 15:28
  • 1
    $\begingroup$ I should also note that generating tuples from a length-$k$ list is isomorphic to getting base-$k$ digits, so: mat = Partition[IntegerDigits[Range[0, 15], 2, 4] + 1, 4] $\endgroup$ – J. M.'s torpor May 5 '20 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.