Add another layout manager demo

This demo uses transforms to place icons on a sphere.
This commit is contained in:
Matthias Clasen 2020-08-27 07:16:38 -04:00
parent fd7a9069ce
commit acc2516b3c
12 changed files with 1135 additions and 2 deletions

View File

@ -136,6 +136,16 @@
<file>demochild.h</file>
<file>demochild.c</file>
</gresource>
<gresource prefix="/layoutmanager2">
<file>demo2layout.h</file>
<file>demo2layout.c</file>
<file>demo2widget.h</file>
<file>demo2widget.c</file>
<file>four_point_transform.h</file>
<file>four_point_transform.c</file>
<file>singular_value_decomposition.h</file>
<file>singular_value_decomposition.c</file>
</gresource>
<gresource prefix="/listview_filebrowser">
<file>listview_filebrowser.ui</file>
<file>listview_filebrowser.css</file>
@ -219,6 +229,7 @@
<file>images.c</file>
<file>infobar.c</file>
<file>layoutmanager.c</file>
<file>layoutmanager2.c</file>
<file>links.c</file>
<file>listbox.c</file>
<file>listbox2.c</file>

View File

@ -0,0 +1,204 @@
#include "demo2layout.h"
#include "four_point_transform.h"
struct _Demo2Layout
{
GtkLayoutManager parent_instance;
float position;
float offset;
};
struct _Demo2LayoutClass
{
GtkLayoutManagerClass parent_class;
};
G_DEFINE_TYPE (Demo2Layout, demo2_layout, GTK_TYPE_LAYOUT_MANAGER)
static void
demo2_layout_measure (GtkLayoutManager *layout_manager,
GtkWidget *widget,
GtkOrientation orientation,
int for_size,
int *minimum,
int *natural,
int *minimum_baseline,
int *natural_baseline)
{
GtkWidget *child;
int minimum_size = 0;
int natural_size = 0;
for (child = gtk_widget_get_first_child (widget);
child != NULL;
child = gtk_widget_get_next_sibling (child))
{
int child_min = 0, child_nat = 0;
if (!gtk_widget_should_layout (child))
continue;
gtk_widget_measure (child, orientation, -1,
&child_min, &child_nat,
NULL, NULL);
minimum_size = MAX (minimum_size, child_min);
natural_size = MAX (natural_size, child_nat);
}
*minimum = minimum_size;
*natural = 3 * natural_size;
}
#define RADIANS(angle) ((angle)*M_PI/180.0);
/* Spherical coordinates */
#define SX(r,t,p) ((r) * sin (t) * cos (p))
#define SZ(r,t,p) ((r) * sin (t) * sin (p))
#define SY(r,t,p) ((r) * cos (t))
static double
map_offset (double x)
{
x = fmod (x, 180.0);
if (x < 0.0)
x += 180.0;
return x;
}
static void
demo2_layout_allocate (GtkLayoutManager *layout_manager,
GtkWidget *widget,
int width,
int height,
int baseline)
{
GtkWidget *child;
GtkRequisition child_req;
int i, j, k;
float x0, y0;
float w, h;
graphene_point3d_t p1, p2, p3, p4;
graphene_point3d_t q1, q2, q3, q4;
double t_1, t_2, p_1, p_2;
double r;
graphene_matrix_t m;
GskTransform *transform;
double position = DEMO2_LAYOUT (layout_manager)->position;
double offset = DEMO2_LAYOUT (layout_manager)->offset;
/* for simplicity, assume all children are the same size */
gtk_widget_get_preferred_size (gtk_widget_get_first_child (widget), &child_req, NULL);
w = child_req.width;
h = child_req.height;
r = 300;
x0 = y0 = 300;
for (child = gtk_widget_get_first_child (widget), i = 0;
child != NULL;
child = gtk_widget_get_next_sibling (child), i++)
{
j = i / 36;
k = i % 36;
gtk_widget_set_child_visible (child, FALSE);
graphene_point3d_init (&p1, 0., 0., 1.);
graphene_point3d_init (&p2, w, 0., 1.);
graphene_point3d_init (&p3, 0., h, 1.);
graphene_point3d_init (&p4, w, h, 1.);
t_1 = RADIANS (map_offset (offset + 10 * j));
t_2 = RADIANS (map_offset (offset + 10 * (j + 1)));
p_1 = RADIANS (position + 10 * k);
p_2 = RADIANS (position + 10 * (k + 1));
if (t_2 < t_1)
continue;
if (SZ (r, t_1, p_1) > 0 ||
SZ (r, t_2, p_1) > 0 ||
SZ (r, t_1, p_2) > 0 ||
SZ (r, t_2, p_2) > 0)
continue;
gtk_widget_set_child_visible (child, TRUE);
graphene_point3d_init (&q1, x0 + SX (r, t_1, p_1), y0 + SY (r, t_1, p_1), SZ (r, t_1, p_1));
graphene_point3d_init (&q2, x0 + SX (r, t_2, p_1), y0 + SY (r, t_2, p_1), SZ (r, t_2, p_1));
graphene_point3d_init (&q3, x0 + SX (r, t_1, p_2), y0 + SY (r, t_1, p_2), SZ (r, t_1, p_2));
graphene_point3d_init (&q4, x0 + SX (r, t_2, p_2), y0 + SY (r, t_2, p_2), SZ (r, t_2, p_2));
/* Get a matrix that moves p1 -> q1, p2 -> q2, ... */
perspective_3d (&p1, &p2, &p3, &p4,
&q1, &q2, &q3, &q4,
&m);
transform = gsk_transform_matrix (NULL, &m);
/* Since our matrix was built for transforming points with z = 1,
* prepend a translation to the z = 1 plane.
*/
transform = gsk_transform_translate_3d (transform,
&GRAPHENE_POINT3D_INIT (0, 0, 1));
gtk_widget_allocate (child, w, h, -1, transform);
}
}
static GtkSizeRequestMode
demo2_layout_get_request_mode (GtkLayoutManager *layout_manager,
GtkWidget *widget)
{
return GTK_SIZE_REQUEST_CONSTANT_SIZE;
}
static void
demo2_layout_class_init (Demo2LayoutClass *klass)
{
GtkLayoutManagerClass *layout_class = GTK_LAYOUT_MANAGER_CLASS (klass);
layout_class->get_request_mode = demo2_layout_get_request_mode;
layout_class->measure = demo2_layout_measure;
layout_class->allocate = demo2_layout_allocate;
}
static void
demo2_layout_init (Demo2Layout *self)
{
}
GtkLayoutManager *
demo2_layout_new (void)
{
return g_object_new (DEMO2_TYPE_LAYOUT, NULL);
}
void
demo2_layout_set_position (Demo2Layout *layout,
float position)
{
layout->position = position;
}
float
demo2_layout_get_position (Demo2Layout *layout)
{
return layout->position;
}
void
demo2_layout_set_offset (Demo2Layout *layout,
float offset)
{
layout->offset = offset;
}
float
demo2_layout_get_offset (Demo2Layout *layout)
{
return layout->offset;
}

View File

@ -0,0 +1,16 @@
#pragma once
#include <gtk/gtk.h>
#define DEMO2_TYPE_LAYOUT (demo2_layout_get_type ())
G_DECLARE_FINAL_TYPE (Demo2Layout, demo2_layout, DEMO2, LAYOUT, GtkLayoutManager)
GtkLayoutManager * demo2_layout_new (void);
void demo2_layout_set_position (Demo2Layout *layout,
float position);
float demo2_layout_get_position (Demo2Layout *layout);
void demo2_layout_set_offset (Demo2Layout *layout,
float offset);
float demo2_layout_get_offset (Demo2Layout *layout);

View File

@ -0,0 +1,172 @@
#include "demo2widget.h"
#include "demo2layout.h"
struct _Demo2Widget
{
GtkWidget parent_instance;
gint64 start_time;
gint64 end_time;
float start_position;
float end_position;
float start_offset;
float end_offset;
gboolean animating;
};
struct _Demo2WidgetClass
{
GtkWidgetClass parent_class;
};
G_DEFINE_TYPE (Demo2Widget, demo2_widget, GTK_TYPE_WIDGET)
static void
demo2_widget_init (Demo2Widget *self)
{
gtk_widget_set_focusable (GTK_WIDGET (self), TRUE);
}
static void
demo2_widget_dispose (GObject *object)
{
GtkWidget *child;
while ((child = gtk_widget_get_first_child (GTK_WIDGET (object))))
gtk_widget_unparent (child);
G_OBJECT_CLASS (demo2_widget_parent_class)->dispose (object);
}
/* From clutter-easing.c, based on Robert Penner's
* infamous easing equations, MIT license.
*/
static double
ease_out_cubic (double t)
{
double p = t - 1;
return p * p * p + 1;
}
static gboolean
update_position (GtkWidget *widget,
GdkFrameClock *clock,
gpointer data)
{
Demo2Widget *self = DEMO2_WIDGET (widget);
Demo2Layout *layout = DEMO2_LAYOUT (gtk_widget_get_layout_manager (widget));
gint64 now;
double t;
now = gdk_frame_clock_get_frame_time (clock);
if (now >= self->end_time)
{
self->animating = FALSE;
return G_SOURCE_REMOVE;
}
t = (now - self->start_time) / (double) (self->end_time - self->start_time);
t = ease_out_cubic (t);
demo2_layout_set_position (layout, self->start_position + t * (self->end_position - self->start_position));
demo2_layout_set_offset (layout, self->start_offset + t * (self->end_offset - self->start_offset));
gtk_widget_queue_allocate (widget);
return G_SOURCE_CONTINUE;
}
static void
rotate_sphere (GtkWidget *widget,
const char *action,
GVariant *parameters)
{
Demo2Widget *self = DEMO2_WIDGET (widget);
Demo2Layout *layout = DEMO2_LAYOUT (gtk_widget_get_layout_manager (widget));
GtkOrientation orientation;
int direction;
g_variant_get (parameters, "(ii)", &orientation, &direction);
self->end_position = self->start_position = demo2_layout_get_position (layout);
self->end_offset = self->start_offset = demo2_layout_get_offset (layout);
if (orientation == GTK_ORIENTATION_HORIZONTAL)
self->end_position += 10 * direction;
else
self->end_offset += 10 * direction;
self->start_time = g_get_monotonic_time ();
self->end_time = self->start_time + 0.5 * G_TIME_SPAN_SECOND;
if (!self->animating)
{
gtk_widget_add_tick_callback (widget, update_position, NULL, NULL);
self->animating = TRUE;
}
}
static void
demo2_widget_snapshot (GtkWidget *widget,
GtkSnapshot *snapshot)
{
GtkWidget *child;
for (child = gtk_widget_get_first_child (widget);
child != NULL;
child = gtk_widget_get_next_sibling (child))
{
/* our layout manager sets this for children that are out of view */
if (!gtk_widget_get_child_visible (child))
continue;
gtk_widget_snapshot_child (widget, child, snapshot);
}
}
static void
demo2_widget_class_init (Demo2WidgetClass *class)
{
GObjectClass *object_class = G_OBJECT_CLASS (class);
GtkWidgetClass *widget_class = GTK_WIDGET_CLASS (class);
object_class->dispose = demo2_widget_dispose;
widget_class->snapshot = demo2_widget_snapshot;
gtk_widget_class_install_action (widget_class, "rotate", "(ii)", rotate_sphere);
gtk_widget_class_add_binding_action (widget_class,
GDK_KEY_Left, 0,
"rotate",
"(ii)", GTK_ORIENTATION_HORIZONTAL, -1);
gtk_widget_class_add_binding_action (widget_class,
GDK_KEY_Right, 0,
"rotate",
"(ii)", GTK_ORIENTATION_HORIZONTAL, 1);
gtk_widget_class_add_binding_action (widget_class,
GDK_KEY_Up, 0,
"rotate",
"(ii)", GTK_ORIENTATION_VERTICAL, 1);
gtk_widget_class_add_binding_action (widget_class,
GDK_KEY_Down, 0,
"rotate",
"(ii)", GTK_ORIENTATION_VERTICAL, -1);
/* here is where we use our custom layout manager */
gtk_widget_class_set_layout_manager_type (widget_class, DEMO2_TYPE_LAYOUT);
}
GtkWidget *
demo2_widget_new (void)
{
return g_object_new (DEMO2_TYPE_WIDGET, NULL);
}
void
demo2_widget_add_child (Demo2Widget *self,
GtkWidget *child)
{
gtk_widget_set_parent (child, GTK_WIDGET (self));
}

View File

@ -0,0 +1,11 @@
#pragma once
#include <gtk/gtk.h>
#define DEMO2_TYPE_WIDGET (demo2_widget_get_type ())
G_DECLARE_FINAL_TYPE (Demo2Widget, demo2_widget, DEMO2, WIDGET, GtkWidget)
GtkWidget * demo2_widget_new (void);
void demo2_widget_add_child (Demo2Widget *self,
GtkWidget *child);

View File

@ -0,0 +1,90 @@
#include "four_point_transform.h"
#include "singular_value_decomposition.h"
/* Make a 4x4 matrix that maps
* e1 -> p1
* e2 -> p3
* e3 -> p3
* (1,1,1,0) -> p4
*/
static void
unit_to (graphene_point3d_t *p1,
graphene_point3d_t *p2,
graphene_point3d_t *p3,
graphene_point3d_t *p4,
graphene_matrix_t *m)
{
graphene_vec3_t v1, v2, v3, v4;
graphene_vec4_t vv1, vv2, vv3, vv4, p;
graphene_matrix_t u, s;
float v[16] = { 0., };
double A[16];
double U[16];
double S[4];
double V[16];
double B[4];
double x[4];
int i, j;
graphene_point3d_to_vec3 (p1, &v1);
graphene_point3d_to_vec3 (p2, &v2);
graphene_point3d_to_vec3 (p3, &v3);
graphene_point3d_to_vec3 (p4, &v4);
graphene_vec4_init_from_vec3 (&vv1, &v1, 1.);
graphene_vec4_init_from_vec3 (&vv2, &v2, 1.);
graphene_vec4_init_from_vec3 (&vv3, &v3, 1.);
graphene_vec4_init_from_vec3 (&vv4, &v4, 1.);
graphene_vec4_init (&p, 0., 0., 0., 1.);
graphene_matrix_init_from_vec4 (&u, &vv1, &vv2, &vv3, &p);
/* solve x * u = vv4 */
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
A[j * 4 + i] = graphene_matrix_get_value (&u, i, j);
B[0] = graphene_vec4_get_x (&vv4);
B[1] = graphene_vec4_get_y (&vv4);
B[2] = graphene_vec4_get_z (&vv4);
B[3] = graphene_vec4_get_w (&vv4);
singular_value_decomposition (A, 4, 4, U, S, V);
singular_value_decomposition_solve (U, S, V, 4, 4, B, x);
v[ 0] = x[0];
v[ 5] = x[1];
v[10] = x[2];
v[15] = 1;
graphene_matrix_init_from_float (&s, (const float *)&v);
graphene_matrix_multiply (&s, &u, m);
}
/* Make a 4x4 matrix that maps
* p1 -> q1
* p2 -> q2
* p3 -> q3
* p4 -> q4
*/
void
perspective_3d (graphene_point3d_t *p1,
graphene_point3d_t *p2,
graphene_point3d_t *p3,
graphene_point3d_t *p4,
graphene_point3d_t *q1,
graphene_point3d_t *q2,
graphene_point3d_t *q3,
graphene_point3d_t *q4,
graphene_matrix_t *m)
{
graphene_matrix_t a, a_inv, b;
unit_to (p1, p2, p3, p4, &a);
unit_to (q1, q2, q3, q4, &b);
graphene_matrix_inverse (&a, &a_inv);
graphene_matrix_multiply (&a_inv, &b, m);
}

View File

@ -0,0 +1,13 @@
#pragma once
#include <graphene.h>
void perspective_3d (graphene_point3d_t *p1,
graphene_point3d_t *p2,
graphene_point3d_t *p3,
graphene_point3d_t *p4,
graphene_point3d_t *q1,
graphene_point3d_t *q2,
graphene_point3d_t *q3,
graphene_point3d_t *q4,
graphene_matrix_t *m);

View File

@ -1,4 +1,4 @@
/* Layout Manager
/* Layout Manager/Transition
*
* This demo shows a simple example of a custom layout manager
* and a widget using it. The layout manager places the children

View File

@ -0,0 +1,198 @@
/* Layout Manager/Transformation
*
* This demo shows how to use transforms in a nontrivial
* way with a custom layout manager. The layout manager places
* icons on a sphere that can be rotated using arrow keys.
*/
#include <gtk/gtk.h>
#include "demo2widget.h"
#include "demo2layout.h"
#include "demochild.h"
GtkWidget *
do_layoutmanager2 (GtkWidget *parent)
{
static GtkWidget *window = NULL;
if (!window)
{
GtkWidget *widget;
GtkWidget *child;
const char *name[] = {
"action-unavailable-symbolic",
"address-book-new-symbolic",
"application-exit-symbolic",
"appointment-new-symbolic",
"bookmark-new-symbolic",
"call-start-symbolic",
"call-stop-symbolic",
"camera-switch-symbolic",
"chat-message-new-symbolic",
"color-select-symbolic",
"contact-new-symbolic",
"document-edit-symbolic",
"document-new-symbolic",
"document-open-recent-symbolic",
"document-open-symbolic",
"document-page-setup-symbolic",
"document-print-preview-symbolic",
"document-print-symbolic",
"document-properties-symbolic",
"document-revert-symbolic-rtl",
"document-revert-symbolic",
"document-save-as-symbolic",
"document-save-symbolic",
"document-send-symbolic",
"edit-clear-all-symbolic",
"edit-clear-symbolic-rtl",
"edit-clear-symbolic",
"edit-copy-symbolic",
"edit-cut-symbolic",
"edit-delete-symbolic",
"edit-find-replace-symbolic",
"edit-find-symbolic",
"edit-paste-symbolic",
"edit-redo-symbolic-rtl",
"edit-redo-symbolic",
"edit-select-all-symbolic",
"edit-select-symbolic",
"edit-undo-symbolic-rtl",
"edit-undo-symbolic",
"error-correct-symbolic",
"find-location-symbolic",
"folder-new-symbolic",
"font-select-symbolic",
"format-indent-less-symbolic-rtl",
"format-indent-less-symbolic",
"format-indent-more-symbolic-rtl",
"format-indent-more-symbolic",
"format-justify-center-symbolic",
"format-justify-fill-symbolic",
"format-justify-left-symbolic",
"format-justify-right-symbolic",
"format-text-bold-symbolic",
"format-text-direction-symbolic-rtl",
"format-text-direction-symbolic",
"format-text-italic-symbolic",
"format-text-strikethrough-symbolic",
"format-text-underline-symbolic",
"go-bottom-symbolic",
"go-down-symbolic",
"go-first-symbolic-rtl",
"go-first-symbolic",
"go-home-symbolic",
"go-jump-symbolic-rtl",
"go-jump-symbolic",
"go-last-symbolic-rtl",
"go-last-symbolic",
"go-next-symbolic-rtl",
"go-next-symbolic",
"go-previous-symbolic-rtl",
"go-previous-symbolic",
"go-top-symbolic",
"go-up-symbolic",
"help-about-symbolic",
"insert-image-symbolic",
"insert-link-symbolic",
"insert-object-symbolic",
"insert-text-symbolic",
"list-add-symbolic",
"list-remove-all-symbolic",
"list-remove-symbolic",
"mail-forward-symbolic",
"mail-mark-important-symbolic",
"mail-mark-junk-symbolic",
"mail-mark-notjunk-symbolic",
"mail-message-new-symbolic",
"mail-reply-all-symbolic",
"mail-reply-sender-symbolic",
"mail-send-receive-symbolic",
"mail-send-symbolic",
"mark-location-symbolic",
"media-eject-symbolic",
"media-playback-pause-symbolic",
"media-playback-start-symbolic",
"media-playback-stop-symbolic",
"media-record-symbolic"
"media-seek-backward-symbolic",
"media-seek-forward-symbolic",
"media-skip-backward-symbolic",
"media-skip-forward-symbolic",
"media-view-subtitles-symbolic",
"object-flip-horizontal-symbolic",
"object-flip-vertical-symbolic",
"object-rotate-left-symbolic",
"object-rotate-right-symbolic",
"object-select-symbolic",
"open-menu-symbolic",
"process-stop-symbolic",
"send-to-symbolic",
"sidebar-hide-symbolic",
"sidebar-show-symbolic",
"star-new-symbolic",
"system-log-out-symbolic",
"system-reboot-symbolic",
"system-run-symbolic",
"system-search-symbolic",
"system-shutdown-symbolic",
"system-switch-user-symbolic",
"tab-new-symbolic",
"tools-check-spelling-symbolic",
"value-decrease-symbolic",
"value-increase-symbolic",
"view-app-grid-symbolic",
"view-conceal-symbolic",
"view-continuous-symbolic",
"view-dual-symbolic",
"view-fullscreen-symbolic",
"view-grid-symbolic",
"view-list-bullet-symbolic",
"view-list-ordered-symbolic",
"view-list-symbolic",
"view-mirror-symbolic",
"view-more-horizontal-symbolic",
"view-more-symbolic",
"view-paged-symbolic",
"view-pin-symbolic",
"view-refresh-symbolic",
"view-restore-symbolic",
"view-reveal-symbolic",
"view-sort-ascending-symbolic",
"view-sort-descending-symbolic",
"zoom-fit-best-symbolic",
"zoom-in-symbolic",
"zoom-original-symbolic",
"zoom-out-symbolic",
};
int i;
window = gtk_window_new ();
gtk_window_set_title (GTK_WINDOW (window), "Layout Manager—Transformation");
gtk_window_set_default_size (GTK_WINDOW (window), 600, 620);
g_object_add_weak_pointer (G_OBJECT (window), (gpointer *)&window);
widget = demo2_widget_new ();
for (i = 0; i < 18 * 36; i++)
{
child = gtk_image_new_from_icon_name (name[i % G_N_ELEMENTS (name)]);
gtk_widget_set_margin_start (child, 4);
gtk_widget_set_margin_end (child, 4);
gtk_widget_set_margin_top (child, 4);
gtk_widget_set_margin_bottom (child, 4);
demo2_widget_add_child (DEMO2_WIDGET (widget), child);
}
gtk_window_set_child (GTK_WINDOW (window), widget);
}
if (!gtk_widget_get_visible (window))
gtk_widget_show (window);
else
gtk_window_destroy (GTK_WINDOW (window));
return window;
}

View File

@ -39,6 +39,7 @@ demos = files([
'images.c',
'infobar.c',
'layoutmanager.c',
'layoutmanager2.c',
'links.c',
'listbox.c',
'listbox2.c',
@ -103,7 +104,11 @@ extra_demo_sources = files(['main.c',
'demotaggedentry.c',
'demochild.c',
'demolayout.c',
'demowidget.c'])
'demowidget.c',
'demo2layout.c',
'singular_value_decomposition.c',
'four_point_transform.c',
'demo2widget.c'])
if harfbuzz_dep.found() and pangoft_dep.found()
demos += files('font_features.c')

View File

@ -0,0 +1,396 @@
#include <string.h>
#include <float.h>
#include <math.h>
#include <glib.h>
/* See Golub and Reinsch,
* "Handbook for Automatic Computation vol II - Linear Algebra",
* Springer, 1971
*/
#define MAX_ITERATION_COUNT 30
static void
householder_reduction (double *A,
int nrows,
int ncols,
double *U,
double *V,
double *diagonal,
double *superdiagonal)
{
int i, j, k, ip1;
double s, s2, si, scale;
double *pu, *pui, *pv, *pvi;
double half_norm_squared;
memcpy (U, A, sizeof (double) * nrows * ncols);
diagonal[0] = 0.0;
s = 0.0;
scale = 0.0;
for (i = 0, pui = U, ip1 = 1;
i < ncols;
pui += ncols, i++, ip1++)
{
superdiagonal[i] = scale * s;
for (j = i, pu = pui, scale = 0.0;
j < nrows;
j++, pu += ncols)
scale += fabs( *(pu + i) );
if (scale > 0.0)
{
for (j = i, pu = pui, s2 = 0.0; j < nrows; j++, pu += ncols)
{
*(pu + i) /= scale;
s2 += *(pu + i) * *(pu + i);
}
s = *(pui + i) < 0.0 ? sqrt (s2) : -sqrt (s2);
half_norm_squared = *(pui + i) * s - s2;
*(pui + i) -= s;
for (j = ip1; j < ncols; j++)
{
for (k = i, si = 0.0, pu = pui; k < nrows; k++, pu += ncols)
si += *(pu + i) * *(pu + j);
si /= half_norm_squared;
for (k = i, pu = pui; k < nrows; k++, pu += ncols)
*(pu + j) += si * *(pu + i);
}
}
for (j = i, pu = pui; j < nrows; j++, pu += ncols)
*(pu + i) *= scale;
diagonal[i] = s * scale;
s = 0.0;
scale = 0.0;
if (i >= nrows || i == ncols - 1)
continue;
for (j = ip1; j < ncols; j++)
scale += fabs (*(pui + j));
if (scale > 0.0)
{
for (j = ip1, s2 = 0.0; j < ncols; j++)
{
*(pui + j) /= scale;
s2 += *(pui + j) * *(pui + j);
}
s = *(pui + ip1) < 0.0 ? sqrt (s2) : -sqrt (s2);
half_norm_squared = *(pui + ip1) * s - s2;
*(pui + ip1) -= s;
for (k = ip1; k < ncols; k++)
superdiagonal[k] = *(pui + k) / half_norm_squared;
if (i < (nrows - 1))
{
for (j = ip1, pu = pui + ncols; j < nrows; j++, pu += ncols)
{
for (k = ip1, si = 0.0; k < ncols; k++)
si += *(pui + k) * *(pu + k);
for (k = ip1; k < ncols; k++)
*(pu + k) += si * superdiagonal[k];
}
}
for (k = ip1; k < ncols; k++)
*(pui + k) *= scale;
}
}
pui = U + ncols * (ncols - 2);
pvi = V + ncols * (ncols - 1);
*(pvi + ncols - 1) = 1.0;
s = superdiagonal[ncols - 1];
pvi -= ncols;
for (i = ncols - 2, ip1 = ncols - 1;
i >= 0;
i--, pui -= ncols, pvi -= ncols, ip1--)
{
if (s != 0.0)
{
pv = pvi + ncols;
for (j = ip1; j < ncols; j++, pv += ncols)
*(pv + i) = ( *(pui + j) / *(pui + ip1) ) / s;
for (j = ip1; j < ncols; j++)
{
si = 0.0;
for (k = ip1, pv = pvi + ncols; k < ncols; k++, pv += ncols)
si += *(pui + k) * *(pv + j);
for (k = ip1, pv = pvi + ncols; k < ncols; k++, pv += ncols)
*(pv + j) += si * *(pv + i);
}
}
pv = pvi + ncols;
for (j = ip1; j < ncols; j++, pv += ncols)
{
*(pvi + j) = 0.0;
*(pv + i) = 0.0;
}
*(pvi + i) = 1.0;
s = superdiagonal[i];
}
pui = U + ncols * (ncols - 1);
for (i = ncols - 1, ip1 = ncols;
i >= 0;
ip1 = i, i--, pui -= ncols)
{
s = diagonal[i];
for (j = ip1; j < ncols; j++)
*(pui + j) = 0.0;
if (s != 0.0)
{
for (j = ip1; j < ncols; j++)
{
si = 0.0;
pu = pui + ncols;
for (k = ip1; k < nrows; k++, pu += ncols)
si += *(pu + i) * *(pu + j);
si = (si / *(pui + i)) / s;
for (k = i, pu = pui; k < nrows; k++, pu += ncols)
*(pu + j) += si * *(pu + i);
}
for (j = i, pu = pui; j < nrows; j++, pu += ncols)
*(pu + i) /= s;
}
else
for (j = i, pu = pui; j < nrows; j++, pu += ncols)
*(pu + i) = 0.0;
*(pui + i) += 1.0;
}
}
static int
givens_reduction (int nrows,
int ncols,
double *U,
double *V,
double *diagonal,
double *superdiagonal)
{
double epsilon;
double c, s;
double f,g,h;
double x,y,z;
double *pu, *pv;
int i,j,k,m;
int rotation_test;
int iteration_count;
for (i = 0, x = 0.0; i < ncols; i++)
{
y = fabs (diagonal[i]) + fabs (superdiagonal[i]);
if (x < y)
x = y;
}
epsilon = x * DBL_EPSILON;
for (k = ncols - 1; k >= 0; k--)
{
iteration_count = 0;
while (1)
{
rotation_test = 1;
for (m = k; m >= 0; m--)
{
if (fabs (superdiagonal[m]) <= epsilon)
{
rotation_test = 0;
break;
}
if (fabs (diagonal[m-1]) <= epsilon)
break;
}
if (rotation_test)
{
c = 0.0;
s = 1.0;
for (i = m; i <= k; i++)
{
f = s * superdiagonal[i];
superdiagonal[i] *= c;
if (fabs (f) <= epsilon)
break;
g = diagonal[i];
h = sqrt (f*f + g*g);
diagonal[i] = h;
c = g / h;
s = -f / h;
for (j = 0, pu = U; j < nrows; j++, pu += ncols)
{
y = *(pu + m - 1);
z = *(pu + i);
*(pu + m - 1 ) = y * c + z * s;
*(pu + i) = -y * s + z * c;
}
}
}
z = diagonal[k];
if (m == k)
{
if (z < 0.0)
{
diagonal[k] = -z;
for (j = 0, pv = V; j < ncols; j++, pv += ncols)
*(pv + k) = - *(pv + k);
}
break;
}
else
{
if (iteration_count >= MAX_ITERATION_COUNT)
return -1;
iteration_count++;
x = diagonal[m];
y = diagonal[k-1];
g = superdiagonal[k-1];
h = superdiagonal[k];
f = ((y - z) * ( y + z ) + (g - h) * (g + h))/(2.0 * h * y);
g = sqrt (f * f + 1.0);
if (f < 0.0)
g = -g;
f = ((x - z) * (x + z) + h * (y / (f + g) - h)) / x;
c = 1.0;
s = 1.0;
for (i = m + 1; i <= k; i++)
{
g = superdiagonal[i];
y = diagonal[i];
h = s * g;
g *= c;
z = sqrt (f * f + h * h);
superdiagonal[i-1] = z;
c = f / z;
s = h / z;
f = x * c + g * s;
g = -x * s + g * c;
h = y * s;
y *= c;
for (j = 0, pv = V; j < ncols; j++, pv += ncols)
{
x = *(pv + i - 1);
z = *(pv + i);
*(pv + i - 1) = x * c + z * s;
*(pv + i) = -x * s + z * c;
}
z = sqrt (f * f + h * h);
diagonal[i - 1] = z;
if (z != 0.0)
{
c = f / z;
s = h / z;
}
f = c * g + s * y;
x = -s * g + c * y;
for (j = 0, pu = U; j < nrows; j++, pu += ncols)
{
y = *(pu + i - 1);
z = *(pu + i);
*(pu + i - 1) = c * y + s * z;
*(pu + i) = -s * y + c * z;
}
}
superdiagonal[m] = 0.0;
superdiagonal[k] = f;
diagonal[k] = x;
}
}
}
return 0;
}
static void
sort_singular_values (int nrows,
int ncols,
double *S,
double *U,
double *V)
{
int i, j, max_index;
double temp;
double *p1, *p2;
for (i = 0; i < ncols - 1; i++)
{
max_index = i;
for (j = i + 1; j < ncols; j++)
if (S[j] > S[max_index])
max_index = j;
if (max_index == i)
continue;
temp = S[i];
S[i] = S[max_index];
S[max_index] = temp;
p1 = U + max_index;
p2 = U + i;
for (j = 0; j < nrows; j++, p1 += ncols, p2 += ncols)
{
temp = *p1;
*p1 = *p2;
*p2 = temp;
}
p1 = V + max_index;
p2 = V + i;
for (j = 0; j < ncols; j++, p1 += ncols, p2 += ncols)
{
temp = *p1;
*p1 = *p2;
*p2 = temp;
}
}
}
int
singular_value_decomposition (double *A,
int nrows,
int ncols,
double *U,
double *S,
double *V)
{
double *superdiagonal;
superdiagonal = g_alloca (sizeof (double) * ncols);
if (nrows < ncols)
return -1;
householder_reduction (A, nrows, ncols, U, V, S, superdiagonal);
if (givens_reduction (nrows, ncols, U, V, S, superdiagonal) < 0)
return -1;
sort_singular_values (nrows, ncols, S, U, V);
return 0;
}
void
singular_value_decomposition_solve (double *U,
double *S,
double *V,
int nrows,
int ncols,
double *B,
double *x)
{
int i, j, k;
double *pu, *pv;
double d;
double tolerance;
tolerance = DBL_EPSILON * S[0] * (double) ncols;
for ( i = 0, pv = V; i < ncols; i++, pv += ncols)
{
x[i] = 0.0;
for (j = 0; j < ncols; j++)
{
if (S[j] > tolerance)
{
for (k = 0, d = 0.0, pu = U; k < nrows; k++, pu += ncols)
d += *(pu + j) * B[k];
x[i] += d * *(pv + j) / S[j];
}
}
}
}

View File

@ -0,0 +1,17 @@
#pragma once
int singular_value_decomposition (double *A,
int nrows,
int ncols,
double *U,
double *S,
double *V);
void singular_value_decomposition_solve (double *U,
double *S,
double *V,
int nrows,
int ncols,
double *B,
double *x);