casadi_kkt.hpp
1 //
2 // MIT No Attribution
3 //
4 // Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl, KU Leuven.
5 //
6 // Permission is hereby granted, free of charge, to any person obtaining a copy of this
7 // software and associated documentation files (the "Software"), to deal in the Software
8 // without restriction, including without limitation the rights to use, copy, modify,
9 // merge, publish, distribute, sublicense, and/or sell copies of the Software, and to
10 // permit persons to whom the Software is furnished to do so.
11 //
12 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
13 // INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
14 // PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
15 // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
16 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
17 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
18 //
19 
20 
21 // SYMBOL "kkt"
22 template<typename T1>
23 void casadi_kkt(const casadi_int* sp_kkt, T1* nz_kkt,
24  const casadi_int* sp_h, const T1* nz_h,
25  const casadi_int* sp_a, const T1* nz_a,
26  const T1* S, const T1* D, T1* w, casadi_int* iw) {
27  // Local variables
28  casadi_int i, k, j, nx, nz;
29  const casadi_int *h_colind, *h_row, *a_colind, *a_row,
30  *kkt_colind, *kkt_row;
31  // Extract sparsities
32  a_row = (a_colind = sp_a + 2) + (nx = sp_a[1]) + 1;
33  h_row = (h_colind = sp_h + 2) + nx + 1;
34  kkt_row = (kkt_colind = sp_kkt + 2) + (nz = sp_kkt[1]) + 1;
35  // Running indices for each row of A
36  for (i = nx; i < nz; ++i) iw[i - nx] = kkt_colind[i];
37  // Reset w to zero
38  casadi_clear(w, nz);
39  // Loop over columns of [H + D_x; A]
40  for (i=0; i<nx; ++i) {
41  // Copy scaled column of H to w
42  for (k=h_colind[i]; k<h_colind[i+1]; ++k) {
43  j = h_row[k];
44  w[j] = nz_h[k] * S[i] * S[j];
45  }
46  // Copy scaled column of A to w
47  for (k=a_colind[i]; k<a_colind[i+1]; ++k) {
48  j = a_row[k] + nx;
49  w[j] = nz_a[k] * S[i] * S[j];
50  }
51  // Add D_x to diagonal
52  w[i] += D[i];
53  // Copy column to KKT
54  for (k=kkt_colind[i]; k<kkt_colind[i+1]; ++k) {
55  j = kkt_row[k];
56  nz_kkt[k] = w[j];
57  if (j >= nx) nz_kkt[iw[j - nx]++] = w[j];
58  }
59  // Zero out w
60  for (k=h_colind[i]; k<h_colind[i+1]; ++k) w[h_row[k]] = 0;
61  for (k=a_colind[i]; k<a_colind[i+1]; ++k) w[a_row[k] + nx] = 0;
62  }
63  // Copy -D_g to diagonal
64  for (i=nx; i<nz; ++i) {
65  nz_kkt[iw[i - nx]++] = -D[i];
66  }
67 }